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orientation of the antenna. Algorithms were included which compare 
the computed field components to the theoretical incident plane 
wave for each stipulated angle of incidence, in order to determine 
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the computational method was attempted by analyzing the 

perturbation indicated for an ideal radome with relative permittivity 
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computed fields which were minimal for axial incidence, but became 
significant for highly canted incidence. 



m 



t.l 

TABLE OF CONTENTS 

I. INTRODUCTION 1 

II. RADOME DESIGN 6 

A . RADOME DESIGN METHODOLOGY 6 

B. RADOME SHAPING FUNCTIONS 8 

C RADOME WALL CONSTRUCTION 1 2 

D. MATERIALS FOR RADOME CONSTRUCTION 16 

III. ELECTROMAGNETIC FIELD THEORY 21 

A. BACKGROUND 21 

B. GENERATION OF THE CORE REGION FIELDS 24 

C CANTING THE PLANAR ARRAY 37 

IV. PROGRAM EXECUTION AND DATA FLOW 41 

A. THE RADOME MODEL DESIGN PROCESS 41 

B. RADSHP EXECUTION 43 

C RADOME CONSTRUCTION VERIFICATION 49 

D. EMRAD EXECUTION 52 

E CORFLD EXECUTION 62 

V. PROGRAM VALIDATION AND RESULTS 73 

VI CONCLUSIONS 93 

APPENDIX A - EMRAD FIELD THEORY 95 

APPENDIX B - RADSHP PROGRAM CODE 110 

APPENDIX C - CORFLD PROGRAM CODE 130 

APPENDIX D - SOFTWARE SOURCES 164 



IV 



LIST OF REFERENCES 165 

INITIAL DISTRIBUTION LIST 167 



v 



I. INTRODUCTION 



Radome design is a composite discipline evolved from 
aeronautical, electrical, mechanical and materials engineering. The 
word radome was coined during World War II as a composite of two 
words, radar and dome [Ref. l:p. ix] . The term radome originally 

described the family of dome-shaped structures designed and built 
to house and protect the radar antennas on aircraft, but now the 
term is applied to any structure that protects a radio frequency 
antenna system. The structure is designed to provide maximum 
protection to the antenna with minimum interference to the 
antenna's electrical performance. 

Radomes on missiles are presently designed to meet many 
diverse specifications. These specification requirements can be 

divided into two categories, electrical and aero-mechanical. Electrical 
design problems include those that have to do with the radome's 
ability to act as a window, ideally transparent, to the electromagnetic 
energy that must pass through the radome to be used by the 
missile’s antenna system. This energy is of a specific frequency 
range and intensity that is within the operating parameters of the 
antenna system. The magnitude and form of the electromagnetic 
energy will be dependent on range from the source, polarization, and 
incident angle. Aero-mechanical design problems address not only 
the design of the structure of the radome, but also the ability of the 
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radome to withstand environmental conditions and battle damage. 
Environmental conditions vary with scenario, and include: normal 
pressure force and viscosity, wind loads, inertial loads, temperature 
variations, moisture resistance, resistance to rain erosion, and other 
diverse factors. Design of a radome for use in battle may require 
that the aero-mechanical design consider resistance to nuclear 
radiation and laser weapons. 

As in any engineering design, there are tradeoffs between the 
various requirements, such that each interrelated component design 
factor must be weighted and considered in the overall system design. 
Numerical models of environmental conditions assist in the analysis 
of the tradeoffs for the aero-mechanical design. Numerical models to 
aid in the electrical design analysis of a radome to protect missile 
antenna systems which operate at high frequencies are available, 
such as the ray-tracing algorithms described by Hirsch and Grove in 
Ref. 2. When a missile antenna system requires operation in the 
resonance region of frequencies for the radome parameters, a 
numerical model is required that can solve the radome interaction 
problem. 

The subject of this thesis is the development of a generalized, 
interactive set of programs to be used in the analysis of the 
interaction between an actual or proposed radome structure and 
incident electromagnetic fields at the lower frequencies of interest, 
where strong resonance characteristics exist. In order to analyze a 
specified radome, a numerical model is constructed using the 
interactive program, RADSHP. This radome generation program 
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provides geometric shaping functions to allow the user to construct a 
multilayer radome model of specified materials. Chapter II provides 
a discussion of radome design methodology, the aerodynamic shaping 
functions available in this package, radome wall constructions, and 
material specification. 

This thesis work is an application of the programs developed 
for use on a personal computer by Morgan [Ref. 3] and Connolly 
[Ref. 4] to calculate electromagnetic scattering by inhomogeneous 
dielectric and/or magnetic bodies of revolution. As discussed in 

Chapter II, radomes are manufactured using dielectric materials in 
order to provide a radio frequency window to the antenna housed in 
the radome's interior. An axisymmetric missile radome structure 

meets the physical structure requirement for rotational symmetry. 
The numerical analysis of the interaction of the specified radome 

structure and the incident electromagnetic energy is performed by 
the code developed by Morgan, EMCAD. Modifications to this code 
resulted in the program EMRAD, which provides data file generation 
to store the modal expansion coefficients which are used by the 
program CORFLD to determine the electromagnetic fields within the 
radome interior. In Chapter III, the theory required for the 

formulation of the computed interior fields of the radome by the 
program CORFLD is discussed. 

The programs developed by Morgan use an optimized 
variational finite-element algorithm in conjunction with a 
tri-regional unimoment method to provide scattering solutions for 

each of multiple incident fields impinging on the specified structure. 
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The programs represent a significant improvement in computational 
speed and versatility over computer codes which find the solution to 
the electromagnetic resonance behavior by use of surface current 
integral equation formulations. An overview of the finite element 
scattering algorithm, along with the concepts embodied in the 
coupled azimuthal potentials and the unimoment method are 
included in Appendix A. Further information is available in 
references 3, 4, 5, 6, and 7. 

Chapter IV discusses the flow of data between the programs 
which have been developed. The general use of the RADSHP 
program to construct the radome data files is discussed, as well as 
the data files generated and their uses. The flow of information into 
and out of the radome analysis program EMRAD is analyzed. The 
execution and operation of CORFLD is discussed. 

In Chapter V, the validation of the interior field computation is 
discussed, wherein the numerically generated fields are compared to 
the incident fields for a radome constructed of air. The results 
produced for radomes constructed of materials presently used in the 
radome industry are presented and discussed. 

The conclusions in Chapter VI will summarize the work 
previously presented, as well as provide recommendations for future 
modifications to the interactive package of programs being used. 

Additional Appendices are included to provide copies of the 
source code for the computer programs, RADSHP and CORFLD, which 
were written as part of this thesis effort. Source listings for 
additional programs, such as EMCADIN, EMESH, and EMRAD are 
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available upon request from Prof. Michael A. Morgan. His address is 
given in Appendix D, Software Sources. 
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II. RADOME DESIGN 



A. RADOME DESIGN METHODOLOGY 

The configuration of a radome for high-speed missiles must 
combine high aerodynamic performance with a highly accurate 
transmission of electromagnetic energy through the radome to the 
missile seeker. As described in Ref. 2, the methodology for 
establishing radome performance requirements is an iterative 
process. Figure 2.1 shows the process as a set of design procedure 
loops that begin with a statement of desired performance, and 
include analytical studies, hardware development, and hardware 
evaluation. These result in a set of optimized radome requirements, 
which is the blueprint from which a radome can be constructed. 
Computer software tools help to maximize the benefit of the work 
done in analytical studies, thereby minimizing the number of 
iterations required for the hardware development and evaluation. 
The programs described in this thesis are designed to analyze the 
interaction between the aerodynamic considerations in radome 
design and the electrical performance of the radome. 

The aerodynamic considerations in the design of a radome 
structure are three-fold: shape-generated drag, thermal stress from 
aerodynamic heating, and structural performance. These three 
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Figure 2. 
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parameters are not independent of each other, and the electrical 
performance of the radome is a strong function of these aerodynamic 
factors. The outer shape of the radome has a large impact on how 
the missile will fly through the air. It is this outer conformation that 
will, along with speed of travel and angle of attack, determine the 
forces that act on the radome structure during flight. The forces 
generated during flight determine the design requirements on the 
radome structure for material strength and thermal stress. 
Fortunately, the shapes that result in the best drag performance are 
those that can also be used to develop radomes that match the 
structural requirements. Unfortunately, the shapes that best meet 
the drag requirements are not necessarily the best for maximizing 
the electrical performance. The programs developed in this thesis 
include the computer aided design (CAD) functions to analyze several 
shaping functions that aerodynamicists have developed to minimize 
the aerodynamic drag forces with respect to radome structure 
parameters. 

B. RADOME SHAPING FUNCTIONS 

An ogive shaping function is one of the most common used in 
radome designs. The ogive shape is popular because it is easily 
fabricated and has relatively acceptable aerodynamic performance at 
high speeds. It is defined [Ref. l:p. 48] as a segment of the arc of a 
circle, whose radius, Fj, is greater than the radius of the radome 

base. A pointed ogive is generated when the circle's radius, Fj, is 
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large enough so that the arc reaches the radome's axis of symmetry. 
The ogive shaping function is described by 



x 




( 2 . 1 ) 



where F 2 and F 3 determine the position of the circular shaping arc, as 
illustrated in Fig. 2.2. The ogive shaping function degenerates to a 
hemisphere when the shaping factors, F 2 and F 3 , are defined as 

F 2 = Fj and F 3 = 0. 

A tangent ogive shaping function generates a radome shape 
such that the surface of the radome is tangent at the base of the 
radome, providing a smooth transition to the missile body [Ref. 8]. A 
tangent ogive shaping function is determined such that 



A shaping function that provides maximum volume-to-drag 
ratio at high flight speeds is the von Karman function [Ref. l:p. 48 j , 




( 2 . 2 ) 



where 




(2.3) 



and 




(2.4) 
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Figure 2.2 Radome shape - nomenclature 



1 0 



where 



F z = cos 1 1 



- (H-z) 
H 



( 2 . 6 ) 



Radomes described by this von Karman shaping function have a 
limited application due to fabrication difficulties related to the 
complex shape. 

In order to minimize drag for a given height-to-radius ratio, a 
Newtonian contour is used. This contour is approximated by the 
power series radome shaping function [Ref. l:p. 49], 



x = R 



f (H-z) ^ F6 
H 



(2.7) 



The exponent, F 6 , provides for shaping. When F 6 = 1.0, the radome 



shape generated is conical. Aerodynamicists normally use F 6 = 0.75 

to approximate the Newtonian contour. 

A shaping function analyzed by Hirsch and Grove [Ref. 2] uses a 
parabolic basis function for the radome shape, 

x 2 = F 7 -z . (2.8) 

The rounded effect developed at the tip of the radome by the 
parabolic shape function can be eliminated by the creation of a 
mathematical "deadzone" for points within a region about the 
radome’s apex. The parabolic shaping coefficient, F 7 , determines the 

extent of the "dead zone" and thereby controls the slope of the 
parabolic surface. The radome shape is determined in terms of the 

shaping coefficient, then scaled to the radome size parameters. The 
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unsealed shaping function, 



x = JF 7 -z , (2.9) 

is used to generate a parabolic radome shape with a height of F 7 . In 
order to eliminate the parabolic rounding of the apex, the "deadzone" 
is established for any points for the unsealed z value within 1.0 of F 7 . 

Within this deadzone, z > ( F 7 - 1 ), the radome shape is determined 
at a line of constant slope, converging at a height of z = F 7 + 1.0. This 
results in the marrying of a pointed apex onto the parabolic radome 
shape, with the magnitude of the shaping coefficient, F 7 , determining 

the relative magnitude of the apex to the total radome height. 

This project contains coding in RADSHP to generate radomes of 
multiple layers for the shaping functions described: the general 

ogive, the tangent ogive, the von Karman shape, the power series 
radome, and the parabolic radome with pointed apex. RADSHP can 
generate a complex radome construction model using different 
shaping functions for each layer. RADSHP generates a material data 
file appropriate for input to the numerical analysis program, EMRAD, 
and a structure data file that can be viewed and modified using 
CURVE DIGITIZER (Appendix 2) , then translated by EMCADIN into an 
input structure data file for EMRAD. 

C. RADOME WALL CONSTRUCTION 

If a material existed which was mechanically strong enough to 

sustain the aerodynamic forces that the radome experiences during 
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flight, and that same material had the electromagnetic properties of a 
vacuum, then the radome could be built to ideally meet all 
requirements. Such a material does not exist, and, in general, as the 
requirement for material strength increases, the electromagnetic 
properties tend to strongly deviate from the properties of the ideal 
vacuum. The transmissivity of the radome tends to decrease as the 
dielectric constants of the construction materials increase. In order 
to meet strength requirements, while minimizing the perturbation 
caused to the electromagnetic field, layer construction of radomes 
has been developed. 

Lamination construction methods allow for an alternating 
combination of thin layers of high density, high strength materials of 
high dielectric constants, with thicker layers of low density, lower 
strength materials of low dielectric constants. Electromagnetic design 
of laminated radomes is more complex than the design of monolithic 
radomes in that there is a greater number of boundary conditions to 
be considered in establishing the lamination layer thicknesses. The 
EMRAD program developed by Morgan [Ref. 3] currently handles 
layered dielectrics with up to five layers. 

Figure 2.3, adapted from the Radome Engineering Handbook 
[Ref. l:p. 46], shows several types of laminated wall construction that 
exhibit the alternation of layers of low dielectric constant and high 
dielectric constant materials. The two-ply sandwich is the simplest 
form of lamination construction with a thin skin of strong, dense 
material of high dielectric constant, supported internally by a thicker 



layer of porous material of low dielectric constant. This design is 
especially attractive at low frequencies as the lateral dimension of 
the thin layer of dense material is only a fraction of the longer 
wavelength. 

The A-sandwich and B-sandwich are two variations of a 
stronger three-layer construction. Although the A-sandwich is 
superior electrically to the B-sandwich, while the placement of the 
high density material in the center of the layer construction will 
strengthen the radome wall by making it stiffer, the porous, low 
density material on the exterior of the radome wall will not sustain 
the aerodynamic forces during flight. 

The B-sandwich's inner and outer skin of high density material 
and the core of low dielectric material provide stiffness and strength 
to the radome structure while providing a dense outer skin. This 
construction has a higher strength-to-weight ratio than a monolithic 
construction. The core thickness is determined by the strength 
requirements and provides spacing between the inner and outer 
skin. Optimal electrical spacing between the skins would by 
frequency dependent at a core thickness of a quarter wavelength 
[Ref l:p. 45]. 

The C-sandwich and multilayer radome wall constructions 
provide for very high strength levels which may be required at high 
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operating speeds or for use in inclement weather. The layer material 
may consist of plastic laminates and ceramics. Multilayered 
construction is attractive where extreme structural rigidity or 
broadband capabilities are primary requirements. [Ref. l:p. 46] 

The radome interacts with the electromagnetic field as a lens, 
modifying the field as it passes through the radome structure. At 
higher frequencies, this interaction can predicted by the use of 
ray-tracing techniques, similar to optical lens analysis [Refs. 2 and 
8]. Studies have been conducted at higher frequencies to design a 
layer of dielectric material within the radome’s interior structure to 
act as a correcting lens [Ref. 8] as shown in Fig. 2.4 . This package of 
programs will allow similar analysis of lens construction at lower 
frequencies, where the radome and lens parameters are within the 
resonance region of the electromagnetic fields. 

D. MATERIALS FOR RADOME CONSTRUCTION 

The choice of material for construction of a radome is another 
engineering design decision. As material science engineers continue 
to develop new materials, particularly low-weight, high-strength 
ceramics, the range of materials available to the radome designer 
increases. The general properties of the three groups of materials 
which are in common use are listed in Ref. 9, and is used here as a 
source for material selection. The groups of materials discussed are 
fiber-reinforced organic resins, refractory oxides, and glass ceramics. 
Physical properties are listed in Table 2.1. 




Figure 2.4 



Interior radome lens of dielectric material 



TABLE 2.1 



RADOME MATERIAL PROPERTIES 
FIBER 

REINFORCED GLASS REFRACTORY 

RESINS CERAMICS OXIDES 



Relative 

permittivity 4-6 5.5 8-9 

(Cr) ' 



Loss tangent 0.01 0.003 0.002 

( tan 8 



Density 



Ultimate 
Tensile 
Strength 
(UTS) 

Compressive 

Strength 20-40 kpsi 35 kpsi 280 kpsi 



Maximum 

Temperature 260 °C 500 °C 1700 °C 



N/A 



160 lb/ft 3 220 lb/ft 3 



40 kpsi 



7 kpsi 



27 kpsi 



In fiber-reinforced organic resins, a variety of resin types are 
used to impregnate fiberglass. Organic resins, such as phenolic, 
polyester, and epoxy are used with the fiber to construct the radome. 
They are then thermally treated to set the structure. Several 

manufacturing techniques are discussed in Ref. 1. Control of layer 
thickness is difficult to maintain and requires finishing to avoid 
aberration. Radomes constructed of these materials have poor 

resistance to rain erosion, but are relatively cheap and easy to 
manufacture, even in large sizes. Maximum temperatures that the 
material can sustain limit use of these materials to missiles with 
speed of less than Mach 1.5 to Mach 2.0. 

Typical glass ceramics used in radome construction are 
Pyrocream and glass-bonded mica. Compared to organic radome 
materials, ceramic materials have much higher resistance to rain 
erosion at supersonic speeds. Ceramics are brittle materials and the 
controls needed during manufacturing to ensure structural integrity 
of the radomes makes them expensive to build [Ref. l:p. 242], 

The refractory oxide used in radome construction is usually 
alumina (A1 2 0 3 ), but magnesia (MgO) and fused silica (Si0 2 ) are also 

used. Refractory oxide manufacturing is also expensive due to the 
high temperatures used in forming these materials to the desired 
radome shape. These materials have high dielectric constants, that 
can undergo a change of up to ten percent between room 
temperature and 1700 °C. Some of the methods by which refractory 
oxide radomes are manufactured are slip casting, isostatic pressing. 



and flame spraying [Ref. 1 :pp. 298-313]. Oxides have a high 
resistance to rain erosion and wall thickness can be controlled to 
within 0.002 inches, allowing for extreme control of aberration. 
Although the dielectric constant of these materials is high, they 
provide the highest degree of thermal protection and place the least 
restriction on missile speed. 
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III. ELECTROMAGNETIC HELD THEORY 



A. BACKGROUND 

The analysis of the interaction between a missile's radome and 
the electromagnetic signal that is radiated or received by that 
missile's radar is a problem similar to those undergoing intensive 
investigation in electromagnetic interaction with other physical 
medium. Ongoing study of medium interactions include other 
military weapons and weapons platforms, forested areas, bodies of 
water, rain, and biological tissues [Refs. 3, 5, and 10]. The study of 
electromagnetic wave penetration through missile radomes allows 
for a special approach to the solution of the problem for the common 
case of rotational symmetry in the construction of the radome. In 
the high-frequency regime, the solution may be found as described 
by Hirsch and Grove [Ref. 2], with application of an asymptotic 
approach using the geometrical theory of diffraction (GTD), which 
incorporates ray-tracing. When lower frequencies are considered, 
and the interaction problem is in the resonance region, a complete 
solution of Maxwell's equations is required. The asymptotic GTD 
method yields unsatisfactory results. 

The solution of Maxwell’s equations for the radome used within 
this project is divided into three regions as illustrated in Fig. 3.1. The 
three regions are the homogenous spherical inner radome core 
(designated Region I); the transition region, where a semi-annular 
conformal mesh is generated for field analysis (designated Region II); 
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Figure 3.1 Three spatial radome regions for 
electromagnetic field analysis 
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and an unbounded exterior scattering region (designated Region III). 
The total solutions for the fields in the core, Region I, and the 
scattered fields in the exterior region, Region III, are produced via 
separate spherical harmonic expansions. Using the unimoment 
method (Appendix A), the modes in these expansions are applied as 
boundary conditions on the inner and outer spherical surfaces 
bounding the mesh and the resultant fields for each such mode are 

found using the finite element method. A similar analysis is 

performed for the incident fields on the outer surface. 

Using superposition, the finite element solutions are then 

combined using the same unknown coefficients as in the original 
expressions. The combined finite element solution is equated, in the 
least squares sense, to the original expressions along the inner and 
outer surfaces. Enforcement of this equality allows solution of the 
expansion coefficients for the core and exterior regions. 

The exterior expansion coefficients are used by EMRAD to 
assemble the scattered fields and bistatic radar cross sections (see 
Refs. 3 and 4). The core region spherical harmonic expansion 

coefficients are used to construct the fields that penetrate the 
radome to illuminate the enclosed antenna. The details of how these 
interior region fields are generated, using the core region spherical 
harmonics, is the topic of the next section. An overview of the 
theory used by Morgan in EMRAD is included in Appendix A. 
Further information is available in Refs. 3, 4, 5, 6, and 7. 
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B. GENERATION OF THE CORE REGION FIELDS 

In this section we will consider the construction of the core 
region fields using the expansion coefficients provided by EMRAD. 
These fields are assembled and plotted in various ways over 
user-selected planar surfaces by CORFLD. 

Within the spherical core, the source-free phasor Maxwell's 
curl equations, with suppressed e^ 1 time variation, are 



-VxE=jcopH 


(3.1) 


VxH=jcoeE , 


(3.2) 



where p and e are the respective permeability and permittivity. The 
p and e will usually be those of free-space for the case of the core 
region in a radome. 

In a homogenous, source-free spherical region, such as the 
radome core, it can be shown [Ref. 11] that the vector 

electromagnetic field can be represented by the combination of two 
types of spherical harmonic fields, transverse magnetic (TM) and 
transverse electric (TE). Using the spherical coordinates of Fig. 3.2, 
these TM and TE fields have the respective properties that H r = 0 and 

E r = 0. Furthermore, the TM and TE vector fields can be generated 
from radially directed vector potentials, such that for the TM fields, 

A = A r r (3.3) 

and 
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Figure 3.2 Spherical coordinates with unit vectors 
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F = F r r , 

where the radial components satisfy 



(3.4) 



A \ 



A r 



r j 



A \ 



+k^ 



A, 



j 



= 0 



(3.5) 



and 



V 


+k 2 


V 


l r J 




l r J 



= 0 . 



(3.6) 



These equations are in the form of the complex scalar 
Helmholtz equation, 



V 2 'P +k 2x ¥ = 0 



(3.7) 



In spherical coordinates, the Helmholtz equation becomes [Ref. 11], 



JLJL 

r 2 9r 



i a f ^ ^ i a 2 ^ 



ay 

9r 



r 2 sin0 96 



Q 9y i 
sin0— - 4 



59 J r 2 sin 2 0 9(J) 2 



+k 2 y=0 . (3.8) 



To find a solution to this differential equation, the method of 
separation of variables is used. Letting 

V = R(r) 0(0) <D(<J>) , (3.9) 

the separation of variables procedure yields the system of ordinary 
differential equations: 

d ^ 



d r 



^pj+[k 2 r 2 -n(n+l)|R =0 , 



(3.10) 
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0 = 0 , 



(3.1 1) 



sine dQ{ de, 



1 d f . d@ 

— - — sin0-— + 



, m 

+ n ( n +1) — — 2 ~ 
sin 0 



and 




d 4>' 



(3.12) 



The R equation is related to Bessel's equation, and its solutions are in 
the form of spherical Bessel functions. Spherical Bessel functions, b n , 
are related to ordinary Bessel functions, B n , as indicated. 



The 0 equation is related to Legendre's equation and its solutions are 
the associated Legendre functions. 



where P n m (cos 0) are the associated Legendre functions of the first 
kind and Q n m (cos 0) are the associated Legendre functions of the 

second kind. The associated Legendre functions of the second kind 
tend toward infinity at points along the z-axis and must be excluded 
for solutions which include this region. 




(3.1 3) 



0(0) = L™(cos0) P™(cos0) , QP(cos0) , 



(3.14) 
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The equation has harmonic sinusoidal type solutions 
d>(<J>) = h(m <)>)—» cos(m<l)), sin(m<|>), e ± -’ m<t> 

= c 1 cos(m<t>)+C2sin(m<}>) 

= C3e +jm<t, + c 4 e _jm ^ . (3.15) 

Together, these separated solutions form the product solutions 

to the Helmholtz equation as 

'J'm.n = b n (kr) P n m (cos0) h(m<))) . (3.16) 

Equations of this type form the elementary spherical harmonic wave 
functions. The general solution to the Helmholtz equation is formed 
by summation of linear combinations of these elementary functions, 
one for the TM fields and one for the TE fields. Each of these 

solutions will be of the form 

'F (r,0><]>) = XX b n (kr) P™(cos6) h(m<J>) . (3.17) 

m n 

The summation over all possible m and n is limited by the 

nature of the Legendre polynomial of the first kind in that P n m = 0 

for n < m. The summation is also limited by the restriction that n 
will not equal zero when m is zero. For both n and m equal to 
zero, then the field evaluated is zero, therefore this solution is also 
excluded. Then m is zero and the summation over n will start 

with n equal to one. This leads to the scalar field expansion 
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'F(r,e,<t>)=S S( c mnCOsm<t) +d mn Sinm<)) P™(cos6) b n (kr). (3.18) 

m=0 n=m 
n*0 

As the Helmholtz equation under consideration is given as 



V 2 ¥ +k 2x ¥ = 0, (3.1 9) 

then for the TM case 

W = ^ , (3.20) 

r 

and for the TE case 

V = . (3.2 1) 

The spherical harmonic solutions of the Helmholtz equations for 
these sets of modes are then of the form 



A r 

Fr 



' = = U kr ) b n (k r) ] P™(cos6) h(m(J)) . 



(3.22) 



In this renormalized form of the solution to the scalar Helmholtz 
equation, the Riccati (or Schelkunoff) spherical Bessel functions are 
present. The Riccati spherical Bessel functions relate to the spherical 
Bessel function and normal Bessel functions as follows: 



j„(kr) = [ (kr) b n (kr) ] = B, , (kr) . 

* 2 ("*2) 



(3.2 3) 



Therefore, in the spherical core region, the solutions for A r and F r are 
of the form 
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= X £( C mn cosm<l) +d mn sinm<|)) P™(cos0) J n (kr) . (3.24) 

m=0 n=m 
n*0 

When expressed in terms of exponential Fourier series in <[>, the 
solutions are of the form 



A r 

F r 



= I Ia mn P”(cos6) J n (kr) . (3.25) 

m- — oo n=m 
n*0 

The coefficients within the two summations are determined by the 
solution of the scattering problem provided by EMRAD. 

The spherical coordinate electric and magnetic field 
components can now be found for the TE and TM cases. For the TM 
case, where 



A r 

F, 



Htm “ Vx (A r r ) 



and 



Etm 



t— VxVx( A r r) , 
jcoe 



the field components are 



E 



r TM 



1 

jcoe 



.2 ^ 

—rr Aj. +k 2 A r 

ar 2 " 1 



Ee 



TM 



i i_ a 2 A r 

jcoe r 3ra0 



F 



i i a 2 A r 

jcoe rsin0 a r a<> 



(3.26) 

(3.27) 

(3.2 8) 

(3.2 9) 
(3.3 0) 
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H r = 0 

r TM 


(3.3 1) 


« _ 1 3A, 

tm rsinG 3 4> 


(3.3 2) 


R _ 1 
n HM r ae • 


(3.3 3) 



For the TE case, where 

HTE=-Vx(F,r) (3.34) 

and 

Bte = jijVxVx(F r r) , (3.35) 

the field components are 



m 

•-» 

II 

O 


(3.3 6) 


E, = - 1 3F - 

TE rsinG 3<> 


(3.37) 


E_ _ 1 aF r 

^TE r 3G 


(3.3 8) 


, ( ? 2 \ 

H r = — fyFr+k 2 ^ 
r TE jcap 3r 2 J 


(3.3 9) 


_ i i a 2 if 

te jcoji r 8r8G 


(3.40) 


ii 1 1 3¥r 

^ jcop rsinG 3r 3<J> 


(3.4 1) 



The total fields are the sum of the TE and TM parts 



31 



(3.4 2) 



E r = 

JO)E 



-\2 

—rrA r +k 2 Aj 

9r 2 r v 



» = L_ aF r , 1 1 ^ A r 

rsin6 a<}) jtoe r 9ra0 



1 



a 2 A t 



a, = I jfi. + J_ 

^ r 90 joe rsin0 9 r 9 <> 



H r = J- 



?t 2 

— -yl> +k 2 F r 
9r 2 r 



(3.4 3) 

(3.44) 

(3.45) 



He = 



1 



rsin0 9<J> 



9A r 1 1 9 2 F r 

- + L 



H - 1 3A r + 1 

t 15- + — 



jo)|i r 9r90 
1 9 2 F r 



(3.46) 

(3.47) 



r d0 jco p. rsin0 9r 9<J) 

The solution for the coefficients in the A r and F r expansions is 

performed by the finite element and unimoment solution method 
explained in Appendix A. The finite element method is used by 
EMRAD and applied to Region II between the inner core and outer 
exterior region. In this finite element region there is the 
inhomogeneous axially symmetric structure of the radome. 

EMRAD generates solutions for incident plane waves with 
Poynting vectors at incident angles specified by the program user, 
such that the Poynting vector's direction is specified by the unit 
propagation vector p = sin a x + cos a z. As any plane wave can be 
decomposed into the sum of two orthogonal, linearly polarized plane 
waves, EMRAD evaluates two plane waves for each incident angle, a 
TM plane wave and a TE plane wave, as shown in Fig. 3.3. The TM 
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z 




Figure 3.3 Incident Poynting vector with TE and TM plane waves 
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and TE plane waves are orthogonal to the Poynting vector and to 
each other. The TE plane wave is always in the y direction. The TM 
wave is in the x-z plane and directed as determined by the 
expression cos ax - sin a z. For incident plane waves specified in 

this manner, the solution for the component values of the E fields in 

the core of the radome can be determined from taylored equations 

for the TM and the TE incident plane waves. For the TM incident 

wave, the core E -field component values are determined as: 



For the TE incident wave, the core E -field component values are 
determined as: 




and 





X ^ simplex b n n(n+l)J n (kr) 




pi 1 1 1 M ^ 

^q(xM)=~y~~^ 'CmSin(m ( t))X a n J 



° r m=0 n=m _ 

n*0 







and 
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With the a n and b n coefficients supplied by EMRAD for each pair of 

plane waves, TM and TE, at each user-specified angle, CORFLD 
generates the interior phasor field components using the equations 
above. The total magnitude of the complex vector field at any point 
within the radome's core is expressed as the square root of the sum 
of the squares of the complex absolute values of the field 
components. This total magnitude is computed and displayed for 
each incident wave by CORFLD, using the a n and b n coefficients 

supplied by EMRAD. 

In order to determine the perturbation of the electromagnetic 
field due to the presence of the radome, CORFLD generates the values 
of the incident plane wave, then compares this to the core region 
field components calculated by EMRAD and CORFLD. The formulas 
for the undisturbed incident TM and TE planes waves are 



The magnitude factor, E 0 = 1.0, is set by EMRAD for both TM and TE 

waves. The unit propagation vector, p, includes the direction of 
propagation for the incident waves and is designated as 



E™=E 0 e j koP * r (xcosa-zsina) 



(3.54) 



and 




(3.5 5) 



p = (x sina+z cosa) . 



(3.5 6) 
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This vector is dotted with the position vector, r, for each point on the 
plane of interest, as defined as 

r = (xx + yy+ zz) . (3.57) 

With the dot product evaluated in the exponential, the TM and TE 
waves are specified as 

E™ =E 0 e -iM<sina«cos«) (xcoso _ 2s j nct) (3.5 8) 

and 



E te =E 0 e -jko(xsina+zcosa) ^ (3 5 9) 

The relative magnitude and phase of the computed field in the 
direction of the theoretical field are calculated and displayed. This 
"scaled dot product" is defined as 

Eda. = ^ . (3-60) 

where E is the computed field, ej is a unit vector pointing in the 

same direction as either the TM or TE incident field specified in 
Equations (3.58) and (3.59), and 

Ei =E 0 e - jko(xsina+zcosa) . (3.61) 

The remaining field term that is computed and displayed by CORFLD 
is the error term. The error term is that portion of the computed 
field which is perpendicular to the incident field vector. In 
normalized form this becomes. 



E 



error 



E-(E»e,)ei 

E[ 



(3.6 2) 
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The relative magnitude of the error term is calculated and displayed 
by CORFLD as a measure of electromagnetic disturbance caused by 
the presence of the radome. This built-in theoretical field 
comparison is used in Chapter V as a measure of modelling accuracy 
by comparing the generated theoretical fields to the fields computed 
and assembled by EMRAD and CORFLD for an ideal radome 
constructed of air (e r = 1). 



C. CANTING THE PLANAR ARRAY 

The planar antenna within a radome can be canted in any 

direction allowed by the gimbal system that is used with the 

antenna. In order to evaluate the field components across a canted 
planar antenna within the radome core, the coordinates of each 
point on the antenna must be converted from the local planar 
coordinate system to the global radome core coordinate system. 

CORFLD uses a transformation from local planar to global core 
coordinates, assuming a two-axis gimbal system where the antenna 

is rotated about the original z-axis, as shown in Fig. 3.4, followed by a 
rotation about the new y-axis, as shown in Fig. 3.5. 

The transformation matrix is formed by multiplying two 
matrices, one for each axis rotation [Ref. 12]. The transformation for 
the rotation about the z-axis is provided by the matrix 



x 



y 



z 



-+COSP 

sinp 

OjO 



Of) 

Of) 

1.0 



-sin(5 

cosp 

0.0 
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(3.6 3) 
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Figure 3.4 Rotation about the z-axis 
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Figure 3.5 Rotation about the new y-axis 
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The transformation for the rotation about the y-axis is provided by 



the matrix 
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y’ 
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cosy 

0D 


0.0 

1.0 


siny 

0J0 
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y" 


(3.64) 
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The complete transformation is given by the product 
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(3.6 5) 



This transformation matrix is used by CORFLD to determine the 
global radome coordinates from the local antenna coordinates for a 
gimbaled antenna rotated by p degrees about the z-axis, then y 
degrees about the new y-axis. With this computational tool the fields 
that are received by an antenna within the core of the radome can be 
determine for any orientation of the antenna and for any computed 
incident field. 
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IV. PROGRAM EXECUTION AND DATA FLOW 



A. THE RADOME MODEL DESIGN PROCESS 

The design process developed by this thesis is represented by 
Fig. 4.1. The radome design procedure is an iterative process of 
selecting radome size and material parameters, radome shaping 
functions, and placement of the axis origin in order to allow for an 
optimal mesh formation by EMRAD. 

The flow of data starts when RADSHP is used to design the 
radome. RADSHP generates a material file which is ready for use by 
EMRAD and a structure data file which is in the format required by 
CURVE DIGITIZER. CURVE DIGITIZER is used to visually check the 
radome structure as designed in RADSHP. After verifying the 
structure data file, the program EMCADIN, developed by Connolly 
[Ref. 4], is used to translate the data file into the form required by 
EMRAD. After data translation, the program EMESH, developed by 
Morgan [Ref. 3], is used to check that the mesh to be generated in 
EMRAD is optimally placed. The placement of the origin in RADSHP 
when generating the structure data file can be modified by a 
placement factor. With the material and structure files are prepared, 
EMRAD is executed to generate the interior coefficients of expansion 
for use by CORFLD. A sample annotated execution through the data 
flow follows. 
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Figure 4.1 Radome model design 



process 
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B. RADSHP EXECUTION 

RADSHP is a menu-driven program that accepts data from the 
keyboard. When the executable version of RADSHP is available, 
RADSHP is invoked by typing 

> RADSHP 



at the DOS prompt. The program will load itself, then provide 
introductory information. 

WELCOME TO RADSHP 

THIS PROGRAM GENERATES RADOME MATERIAL AND SHAPE DATA 
FILES FOR USE WITH EMRAD, BASED ON USER INPUT. 

RADSHP prompts for data file names for the material and structure 
output files. Sample keyboard inputs are provided. User inputs are 
preceded by "==>". 

TWO FILES ARE REQUIRED TO HOLD PROGRAM OUTPUT DATA. 

PLEASE ENTER THE NAME OF THE FILE TO BE USED TO 
HOLD THE PROGRAM MATERIAL PARAMETER DATA. THE EXTENSION 
.MAT WILL BE ADDED AUTOMATICALLY, I. E. Y OURNAME.MAT 

==> MATFIL 

PLEASE ENTER THE NAME OF THE FILE TO BE USED 
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TO HOLD THE PROGRAM OUTPUT DATA. THE EXTENSION .DAT WILL 
WILL BE ADDED AUTOMATICALLY, I. E. YOURNAME.DAT 



==> STRFIL 

RADSHP will allow the construction of radomes with up to five layers 
of construction. At present, each of these layers must completely 
enclose the previous layered region. 

ENTER THE NUMBER OF LAYERS (.LE. 5) 

==> 2 

RADSHP prompts for the parameters of each layer of construction. 

Each layer of the radome will have material constants, 
an outer height of the radome layer, H, 
a outer base radius, R, and 

a shape function, & shape factor parameters (if appl.), 

The first "layer" will be the inner core of the radome, 
normally filled with air. 

Layer No: 1 

This is the inner core of the radome 



Please enter the radome layer height, H, 
and the outer base radius, R, for this layer. 

==> 1.25,0.50 



44 



RADSHP provides a menu of the shaping functions available for 
radome construction. The curved parabola is not normally used for 
outer layers of the radome, but is included as an option for inner 
layer construction. 

The following shaping functions are available, 

1. General Ogive 

2. Tangent Ogive 

3. von Karman 

4. Power Series 

5. Parabola, curved 

6. Parabola, capped by pointed apex 

Indicate desired shaping function by entering a number 
1,2, ... or 6 



==> 2 



Each shaping function requires input of their respective 
shaping factors. These are prompted for and then listed for each 
layer as the radome is constructed. The tangent ogive and the von 
Karman functions have shaping factors that are determined by the 
height and base radius of the radome layer; these factors are 
determined and listed. The factor designation follows the form 
adopted in the equation in Chapter II. 

The tangent ogive shaping function determines two 
factors for the radome construction. These are calculated 
based on the height and radius of the radome layer. 

These factors are F4 = 6.50 and F5 = 6.38 
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LAYER SHAPE 



HEIGHT RADIUS FI 



1 TAN Ogive 1.25 .50 .00 



F2 


F3 


F4 


F5 


F6 


F7 


.00 


.00 


6.50 


6.38 


.00 


.00 



Layer No: 2 

This is the outer layer of the radome 



Please enter the radome layer height, H, 
and the outer base radius, R, for this layer. 

==> 1.30,0.55 

The following shaping functions are available, 

1 . General Ogive 

2. Tangent Ogive 

3. von Karman 

4. Power Series 

5. Parabola, curved 

6. Parabola, capped by pointed apex 

Indicate desired shaping function by entering a number 
1,2, ... or 6 



The tangent ogive shaping function determines two 
factors for the radome construction. These are calculated 
based on the height and radius of the radome layer. 

These factors are F4 = 6.42 and F5 = 6.28 
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LAYER 


SHAPE 


HEIGHT 


RADIUS 


Fl 


F2 


F3 


F4 


F5 


F6 


F7 


1 


TAN Ogive 


1.25 


.50 


.00 


.00 


.00 


6.50 


6.38 


.00 


.00 


2 


TAN Ogive 


1.30 


.55 


.00 


.00 


.00 


6.42 


6.28 


.00 


O 

q 



After all the structure forming data is entered, RADSHP 
prompts for the material relative permittivity parameter for each 
layer of construction. Layer 1 is the inner-most region and the 
permittivity will normally equal unity. 

Layer No: 1 

This is the inner core of the radome 
Enter Er for this Layer: 

==> 1.0 

Layer No: 2 

This is the outer layer of the radome 
Enter Er for this Layer: 

==> 4.0 

With the structure and material parameters entered, RADSHP 
prompts for headers to identify the radome constructed. This header 
information is stored in the material data file and will follow the data 
files throughout the execution of the programs in this package. 
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THE HEADERS ALLOW THE USER TO IDENTIFY THIS 
SET OF DATA FROM ALL OTHER SETS. 

PLEASE ENTER HEADER#! (64 CHARACTERS MAX) 



==> This is a test data set for RADSHP example 



PLEASE ENTER HEADER #2 (64 CHARACTERS MAX) 
==> Header info stays with data set throughout 



With all the data input and calculations completed, 
reports the generation of the material and structure 
terminates. 

Computations and output are completed, the data file, MA TELL. MAT 
holds the material specifications for input to EMCAD. 

The data file, STRFIL.DAT, holds the radome shape data for 
viewing in the CAD package, CURVE DIGITIZER. 

The structure file required for input to EMCAD can be 
generated by the executing the program EMCADEN, 
using, STRFIL.DAT, as input to EMCADIN, 

The output file from EMCADIN will be in proper format 
for input to EMCAD for the structure file. 

Stop - Program terminated. 



> 



RADSHP 
files and 
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C. RADOME CONSTRUCTION VERIFICATION. 



The construction of the radome by RADSHP is verified by 
entering the CURVE DIGITIZER package and displaying the structure 
data file generated by RADSHP. Figure 4.2 is a copy of the CURVE 
DIGITIZER display of the radome shape created in the example 
execution above. Upon review and acceptance of the structure data 
file, EMCADIN is used to translate the data file into the form required 
by EMRAD. EMCADIN is invoked by typing 

>EMCADIN 

at the system prompt. The program will load itself and displays the 
greeting heading and prompts for the structure data file as input. 



WELCOME TO EMCADIN 

THIS PROGRAM ACCEPTS PROPERLY FORMATTED INPUT DATA 
FROM CURVE DIGITIZER AND CONVERTS IT TO A FORM 
WHICH CAN BE USED BY EMCAD. 

PLEASE PRESS ENTER TO CONTINUE. 

INPUT FILENAME 

INPUTFN IS THE INPUT FILENAME OF THE FILE CONTAINING 
THE DATA OF INTEREST WHICH WILL BE INTERPOLATED ON AND 
CONVERTED TO A FORM WHICH CAN SUBSEQUENTLY BE USED BY 
EMCAD 

PLEASE ENTER THE NAME OF THE INPUT FILE. THE EXTENSION 
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z-axis 




Figure 4.2 Tangent ogive created by RADFLD 
as displayed by CURVE DIGITIZER 
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OF THE FILENAME MUST BE .DAT, I.E. FILENAME.DAT 



==> STRFIL.DAT 



EMCADIN then prompts the user for the data file name to hold the 
structure file that is prepared for EMRAD execution. 

OUTPUT FILENAME 

OUTPUTFN IS THE OUTPUT FILENAME OF THE FILE CONTAINING 

THE DATA WHICH HAS BEEN INTERPOLATED AND CONVERTED TO A 

FORM WHICH WILL BE USED BY EMCAD. 

PLEASE ENTER THE NAME OF THE OUTPUT FILE. THE EXTENSION 

OF THE FILENAME MUST BE INCLUDED I.E. FILENAME.DAT 

==> STRFIL.STR 

There are several data files used within these packages, the 
default file name extensions should be: ".MAT" for material data files, 
".DAT" for structure data files suitable for CURVE DIGITIZER, and 
".STR" for structure files suitable for EMRAD. 

After receiving the data file names to work with, EMCADIN 
prompts for angular resolution requirements. EMCADIN uses linear 
interpolation to provide the capability to generate output data with 
better resolution than the input data. The default resolution of data 
generation by RADSHP is one thousand points in the x-z plane, using 
180 points of angular resolution in EMCADIN is sufficient for good 
results. 
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RESOLUTION 



DELTHE IS THE USER INPUT VALUE OF THE DESIRED THETA 
RESOLUTION IN DEGREES. 

PLEASE ENTER THE DESIRED DELTA THETA VALUE IN DEGREES. 

==> 180. 

Upon completion of the data file translation by EMCADIN, the 
data files required for execution by EMRAD are complete and in the 
proper format. MATFIL.MAT is the example material file and 
STRFIL.STR is the example structure file. Before executing EMRAD, 
the efficiency of mesh generation to be expected during EMRAD can 
be previewed by using the program EMESH developed by Morgan. 
Figure 4.3 shows the mesh generated by EMESH for the example 
radome, as displayed by CURVE DIGITIZER. 

D. EMRAD EXECUTION 

With the data files generated and an optimal mesh solution, the 
data flow is ready to initiate EMRAD. EMRAD is a menu-driven 
program developed by Connolly and Morgan and has been slightly 
modified to provide the interior spherical coefficients of expansion 
for use by CORFLD. EMRAD can be invoked by typing 

>EMCAD 

at the DOS prompt. The program will load itself and provide 
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Figure 4. 3 Tangential ogive with mesh generated by EMESH 
as displayed by CURVE DIGITIZER 
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introductory information. 



****** WELCOME TO EMRAD ****** 



THIS PROGRAM COMPUTES EM FIELD SCATTERING FROM PENETRABLE 
BODIES OF REVOLUTION USING THE COUPLED AZIMUTHAL POTENTIAL 
(CAP) FORMULATION IN CONJUNCTION WITH A VARIATIONAL FINITE- 
ELEMENT TECHNIQUE AND A TRI-REGIONAL UNIMOMENT METHOD. THE 
NUMERICAL SOLUTION IS PERFORMED USING A TWO-SWEEP RICCATI 
TRANSFORM WITH CIRCUMFERENTIAL MARCHING. 

UPDATED VERSION OF ORIGINAL U.C. BERKELEY CODE. 

MODIFICATIONS BY M.A. MORGAN MAR 1987 - DEC 1989 

E.M. CONNOLLY JUL 1987 - APR 1988 
R.J. VINCE SEP 1989 - DEC 1989 . 



PLEASE PRESS ANY KEY TO CONTINUE. 



******* EMRAD INPUT ******* 



EMRAD ALLOWS A USER TO INPUT THE NECESSARY INFORMATION 
ACCORDING TO THE LEVEL OF HIS EXPERTISE. 

1 NOVICE LEVEL - GIVES BRIEF EXPLANATIONS/DESCRIPTIONS 

OF INPUT VALUES. DESCRIBES THE REQUIRED 
FORMAT FOR THE INPUT VALUE. GIVES 
TYPICAL VALUES WHERE APPLICABLE. 

2 EXPERT LEVEL - ASSUMES THE USER IS FAMILIAR WITH THE 

INPUT PARAMETERS AND FORMATS. SIMPLY 
PROMPTS FOR REQUIRED INPUTS. 

3 EXIT EMRAD 
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PLEASE SELECT THE LEVEL OF EXPERTISE BY ENTERING 1, 2, OR 3. 



==> 1 



The names of the data files generated by RADSHP and translated by 
EMCADIN are entered at the keyboard when prompted by EMRAD. 



EMRAD REQUIRES DATA INPUT FROM TWO FILES PREPARED 
IN ACCORDANCE WITH THE USERS MANUAL. THESE FILES MUST 
CONTAIN MATERIAL PARAMETERS AND STRUCTURE PARAMETERS, 
RESPECTIVELY. 

ENTER MATERIAL PARAMETER FILE (DiFTLENAME.EXTENSION) : 



==> MATFIL.MAT 



ENTER STRUCTURE DATA FILE (D:FILENAME.EXTENSION): 

==> STRFIL.STR 

EMRAD opens both files and displays the header information which 
RADSHP placed in the material file. 

This is a test dataset for RADSHP example 
Header will stay with dataset throughout 

EMRAD creates several datafiles. The first will be the diary file for 
the execution of EMRAD. This file will have the ".OUT" extension. 
The second data file generated will have the same root name as the 
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diary file, this is the data file which holds the spherical coefficients of 
expansion for the core field generation. This second file will have the 
extension ".INT" and will be used as input to CORFLD. 



THE OUTPUT DATA FILE IS THE FILE CONTAINING THE OUTPUT 
RESULTS OF ALL NUMERICAL CALCULATIONS CONDUCTED BY EMRAD. 
THE FORMAT FOR THIS INPUT IS FILENAME ONLY. NO EXTENSION IS 
REQUIRED OR DESIRED. EMRAD AUTOMATICALLY APPENDS AN 
EXTENSION OF .OUT TO YOUR FILENAME. 

A SECOND OUTPUT DATA FILE WILL BE CREATED TO STORE THE 
COEFFICIENTS FOR THE INTERIOR CORE EXPANSION. EMRAD WILL 
APPEND AN EXTENSION OF .INT FOR THIS DATA FILE. 

PLEASE ENTER THE FILENAME OF THE OUTPUT DATA FILES. 

==> INTFIL 



EMRAD was developed by Morgan and Connolly to determine the 
scattered field and the radar cross section of the body of revolution 
under study. Four output data files are created for input to separate 
graphing routines. 

THE OUTPUT GRAPHICS DATA FILE IS THE FILE CONTAINING THE 
OUTPUT DATA FROM EMRAD TO BE USED AS INPUT TO GRAPHING 
ROUTINES. THE FORMAT FOR THIS INPUT IS FILENAME ONLY. NO 
EXTENSION IS REQUIRED OR DESIRED. EMRAD AUTOMATICALLY 
APPENDS AN EXTENSION OF .TMT, .TMP, .TET, AND .TEP TO YOUR 
FILENAME AS IT PRODUCES FOUR OUTPUT FILES FOR GRAPHICS. 

.TMT - — > TM INCIDENCE, F-THETA 
.TMP - — > TM INCIDENCE, F-PHI 



56 



.TET - — > TE INCIDENCE, F-THETA 
.TET - — > TE INCIDENCE, F-PHI 



PLEASE ENTER THE FILE NAME OF THE OUTPUT DATA FILE. 

==> SCATGRF 

EMRAD prompts for a graphics caption to be included with the data 
set, this caption is passed along to CORFLD and is included with each 
plotted figure. This caption affords the user the ability to identify 
this data set from all other sets. 

THE GRAPHICS CAPTION IS A PERSONALIZED CAPTION ALLOWING THE 
USER TO IDENTIFY THIS SET OF GRAPHS FROM ALL OTHER SETS. 

THE MAXIMUM LENGTH OF THIS CAPTION IS 64 CHARACTERS. 

NOTE: WHEN USED WITH THE GRAPHICS PACKAGE, THE PROGRAM IS 
ABLE TO DIFFERENTIATE BETWEEN UPPER CASE AND LOWER CASE 
CHARACTERS. 

PLEASE ENTER ANY GRAPHICS CAPTION YOU DESIRE 
==> Graphics caption 

As discussed in Appendix A, EMRAD utilizes a semi-annular 
conformal mesh in order to conduct electromagnetic field 
calculations. EMRAD prompts the user for a minimum and maximum 
mesh density. EMRAD conducts the analysis of the radome structure 
using a linear variation of density from maximum density decreasing 
to minimum density. 

DMIN AND DMAX ARE PARAMETERS OF MESH DENSITY IN TERMS OF 
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ELEMENTS/INTERIOR LAMBDA. INPUT VALUES ARE EXPECTED TO BE 
REAL, I.E. THE DECIMAL POINT MUST BE INCLUDED. 

TYPICAL VALUES ARE DMIN = 10. AND DMAX = 15. 

PLEASE INPUT DMIN 

==> 15 



PLEASE INPUT DMAX 
==> 18 

EMRAD allows for analysis of up to five incident fields. EMRAD 
prompts for the number of incident angles, and then for the value of 
each angle, in degrees. EMRAD uses the poynting vector direction for 
the incident angle. This translates to an angle of 180 degrees for an 
angle incident at the radome’s appex, 135 degrees for an angle 
incident 45 degrees off the apex, and 90 degrees for an incident 
angle off the beam of the radome. 

THE NUMBER OF INCIDENT FIELD ANGLES IS THE TOTAL NUMBER OF 
INCIDENT FIELDS THAT IMPINGE ON THE OBJECT OF INTEREST. 

THIS PROGRAM ALLOWS A MAXIMUM OF FIVE INCIDENT FIELD 
ANGLES. WHEN ENTERING YOUR ANSWER PLEASE DO NOT INCLUDE A 
DECIMAL BECAUSE THE INPUT MUST BE IN INTEGER FORMAT, I.E. 



NA = 3 , NA = 1 



PLEASE INPUT THE NUMBER OF INCIDENT FIELD ANGLES. 
==> 2 
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ENTER INC FLD ANGLE (DEG) FOR # 1 



==> 180 

ENTER INC FLD ANGLE (DEG) FOR # 2 
==> 135 

EMRAD performs scattering calculations by starting at 0=0, and 
progressing clockwise to 0 = 180 in equal increments, the menu 
allows the user to determine the size of the increments. 

THE NUMBER OF SCATTERING FIELD THETA POINTS DETERMINES THE 
SPACING BETWEEN THETA POINTS DURING EMRAD ITERATIONS AND 
CALCULATIONS. 

DELTA THETA = 180/ (NUMBER THETA POINTS - 1) SO... 

NUMBER THETA POINTS = 37 ----- > DELTA THETA = 5 DEGREES 

NUMBER THETA POINTS = 19 > DELTA THETA = 10 DEGREES 

WHEN ENTERING YOUR ANSWER, PLEASE DO NOT INCLUDE A DECIMAL 
BECAUSE THE INPUT MUST BE IN INTEGER FORMAT. I.E. NT = 19 

PLEASE INPUT THE NUMBER OF SCATTERING FIELD THETA POINTS 
==> 19 

As part of the scattering analysis, EMRAD allows for eight angles of 
user observation. 
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THE NUMBER OF PHI ANGLES IS THE TOTAL NUMBER OF PHI ANGLES. 
THIS PROGRAM ALLOWS A MAXIMUM OF THREE PHI ANGLES. WHEN 
ENTERING YOUR ANSWER PLEASE DO NOT INCLUDE A DECIMAL 
BECAUSE THE INPUT MUST BE IN INTEGER FORMAT. I.E. 



NP = 3 , NP = 1 



PLEASE INPUT THE NUMBER OF PHI ANGLES. 

==> 0 

EMRAD truncates the modal analysis of the electromagnetic field 
infinite fourier expansions for the incident plane waves and the 
interior core as well as the exterior scattered fields. Mathematical 
estimates of truncations are provided and the user is prompted for 
truncation limits. Higher truncations may increase the modelling 
accuracy. 

ENTER MSTOP (.LE. 13) 

ESTIMATED "MINIMUM" VALUE IS: 4 

==> 5 



Modal truncations for the number of spherical harmonic modes in 
the spherical radome core and external scattering region are 
estimated. The user is prompted for truncation limits for the 
spherical harmonics. 
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ENTERING NO. OF INTERNAL AND EXTERNAL EXPANSION MODES 



ESTIMATED (KPRMIN) "MINIMUM" INTERNAL N1 = 4 

ENTER N1 (.GE. MSTOP): 



ESTIMATED (KO*RMIN) "MINIMUM" EXTERNAL N2 = 9 

ENTER N2 (.GE. MSTOP): 

==> 9 

Check that N1+N2 is .LE. 35 .... Otherwise Abort 



==> [RETURN] 



When executing the numerical solution, EMRAD requires a large 
amount of temporary disk storage space. The amount of storage 
required is estimated and provided to the user. The user must select 
a disk drive that has adequate memory available for usage. The use 
of a ram drive configuration can increase the speed of execution by 
EMRAD. 

***** Temporary Storage of Riccati Matrices ***** 

Estimated Disk Space Needed in MegaBytes: 1.930000E-KX) 

ENTER Hard Drive or RAM Disk Letter ( C, D, E ... ): 

==> E 
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With all of the data and parameters specified, EMRAD begins the 
numerical analysis of the electromagnetic interaction of the radome 
structure and the specified incident plane waves. Execution is 
monitored by reports to the terminal screen. 

INDEXING PROGRAM THROUGH VALUES OF M 

M-LOOP .... M = 0 

CALL MESH 

ENTER I-LOOP; NO STEPS TO COMPLETE: 38 

1= 1 NO STEPS TO GO: 37 

CALLLODER 
CALL MARCH 

1= 2 NO STEPS TO GO: 36 

EX M-LOOP, DATOUT 

* ******* ******* emrad COMPLETED ************** 

Stop - Program terminated. 

Upon termination of the program EMRAD, the coefficients for analysis 
of the radome core fields have been stored in the data file, ready for 
CORFLD execution. 

E. CORFLD EXECUTION 

CORFLD is a menu driven program that accepts data from the 
data file created by EMRAD and from the keyboard. When the 
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executable version of CORFLD is available, CORFLD is invoked by 
typing 

> CORFLD 

at the DOS prompt. The program will load itself, then provide 
introductory information. 

WELCOME TO CORFLD 

CORFLD ASSEMBLES THE FIELD ON THE INTERIOR OF A RADOME 
WITH THE COEFFICIENTS SUPPLIED BY AN EMRAD OUTPUT FILE. 



CORFLD prompts for the name of the data file generated by EMRAD. 
Only the root file name is required, CORFLD attaches the ".INT” to the 
root file name. 



THE OUTPUT DATA FILE IS THE FILE CONTAINING THE OUTPUT 
RESULTS OF INTERIOR FIELD CALCULATIONS CONDUCTED BY EMRAD. 
THE FORMAT FOR THIS INPUT IS FILENAME ONLY. 

NO EXTENSION IS REQUIRED OR DESIRED. 

CORFLD AUTOMATICALLY APPENDS AN EXTENSION OF .INT TO YOUR 
FILENAME. 

PLEASE ENTER THE FILENAME OF THE OUTPUT DATA FILE. 



==> INTFIL 



Header information that has followed the data flow is obtained from 
the ".INT" data file by CORFLD and displayed on the screen. Notice 
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that the graphics caption entered during EMRAD execution as well as 
the headers entered during RADSHP execution are displayed. 

TEST DATA FOR USER EXAMPLE 

THIS HEADING FOLLOWS THE DATA THROUGHOUT 

Graphics caption 

CORFLD reads the material parameters stored in the ".INT" file and 
displays them on the screen. 



SRE= (1. 000000, -5.000000E-07) 

SRU= (1. 000000, -5.000000E-07) 

KR = (1 .000000,- 1 .000000E-06) 

CORFLD reads the number of incident angles evaluated by EMRAD 
and displays each angle as the poynting vector angle of incidence and 
also as the angle of incidence to the radome exterior. The radome 
incident angle is displayed so that the user can use the angle as a 
reference when canting the planar antenna. 

Incoming Inc. Angle # 1 

EMRAD Poynting angle of incidence = 1 80.000000 

Radome incident angle = 0.000000E+00 

Incoming Inc. Angle # 2 

EMRAD Poynting angle of incidence = 135.000000 

Radome incident angle = 45.000000 

After displaying the angles of incidence available for consideration, 
EMRAD prompts the user for the angle of interest. Only one angle 
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may be selected for each execution of CORFLD. In order to consider 
other angles, CORFLD must be re-executed by the user. 

Select Angle of Interest By Number 
=> 1 

The angle of consideration is angle # 1 = 0.000000E+00 degrees 

After an angle is selected for analysis, the user may enable a rotation 
of the planar array. 

The planar array is set by default to lie in the 
x-y plane, with boresight in the z-direction. 

Do you want to change the array boresight ? 

( "Y" or "N" ) 



==> N 

CORFLD reads the data from the ".INT" file and assembles the interior 
field components. 

MSTART= 0 MSTOP= 5 
Reading Coefficients and Constructing 
Fields for Each Azimuthal "m" Mode 

Terminating M-Loop at m=l for +/- Z-Axis Incidence 
Azimuthal Mode "m" = 0 
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Reading Coefficients 



Azimuthal Mode "m" = 1 



Reading Coefficients 
Generate Fields Across Array 
SPHERICAL FIELDS COMPUTED 



When the interior fields have been computed, the user may check for 
the field component values for any point on the planar array. 



Check Values of Coordinates and E-Field ? (Y/N): 
==> Y 

Enter Array Triple Index to Check (J,L,I) 

J= x-point 1,2,.., 50 

L= y-point 1,2,.., 50 

1= 1 for TM, 2 for TE 



XP,YP,ZP: -5.115592E-02 -5.115592E-02 0.000000E+00 

XC,YC,ZC: -5.1 15592E-02 -5.115592E-02 0.000000E+00 

IE-RADI = 5.023 160E-01 
IE-THI = 1.391270E-02 
IE-PHII = 5.022836E-01 
IE-TOTI = 7.104954E-01 



Look at Another Point ? (Y/N): 
==> N 



The user is supplied with the selection menu to allow access to any 
component representation of the interior fields across the planar 
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antenna which is oriented as specified by the user. Component phase 
and magnitudes of the interior fields as well as the total field 
magnitudes are available, 



CORFLD OUTPUT GRAPHICS PRESENTATION 
CORE FIELD REPRESENTATION SELECTION MENU 



Please select the type of field representation of interest 
CORE will display the following field representations: 



1. TM INCIDENT FIELD, E-THETA 

2. TM INCIDENT FIELD, E-PHI 

3. TM INCIDENT FIELD, E-RADIAL 

4. TM INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE 

5. TE INCIDENT FIELD, E-THETA 

6. TE INCIDENT FIELD, E-PHI 

7. TE INCIDENT FIELD, E-RADIAL 

8. TE INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE 

9. COMPARE COMPUTED FIELD TO INCIDENT FIELD 

10. CHANGE ASPECT RATIO 

0. FINISHED WITH THIS ANGLE OF INCIDENCE 



Indicate your selection entering "1", "2",.. .,"10", or "0" 
==> 4 



When the user has made a selection from the menu, the three 
dimensional (3-D) plotting routine is called and the user is allowed to 
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adjust the view angle for best observation of the 3-D plot, see 
Fig. 4.4, 

CALL PLOTTING ROUTINE 

A black and white 'stick' representation as well as a color 3-D fill 
representation are available to the user. The 'stick' representation 
generates the clearest output for printing and is used for all the data 
presented here. 

Enter 0 for 3D-stick or 1 for 3D-fill option 
==> 0 

(** SEE FIGURE 4.4 **) 

Default view angles are: phi=-45.,theta=70. 

ENTER NEW ANGLES PHI & THETA (DEG) OR (9,9) WHEN DONE 
==> 9,9 

When the user is finished with an electric field component, the new 
view angle "9,9" is entered and the user is returned to the core field 
selection menu. 

CORE FIELD REPRESENTATION SELECTION MENU 

9. COMPARE COMPUTED FIELD TO INCIDENT FIELD 
Indicate your selection entering "1", "2",.. .,"10", or "0" 
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The datafile used to generate fields is intfil.INT 
TEST DATA FOR USER EXAMPLE 
THIS HEADING FOLLOWS THE DATA THROUGHOUT 



Graphics caption 

TM INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE 




Min. Z = +.631; Max. Z = +.825; average Z = +.656 
Inc. A = +.000; Plane G = +.000; Uieu: phi = + 45.0 theta= + 75.0 



Figure 4.4 3-D plot of E-total field magnitude for example ogive radome 
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When the user wishes to see the amount of perturbation caused to 
the incident plane wave by the radome, selection "9" is entered, 



==> 9 



Selection "9" sends the user into Menu 2, for core field comparison to 
the ideal theoretical plane wave which would be incident on the 
antenna without the radome present. 



CORE FIELD REPRESENTATION SELECTION MENU 2 



Please select the type of field representation of interest 
CORE will display the following field representations: 



1. TM INCIDENT FIELD, 

2. TM INCIDENT FIELD, 

3. TE INCIDENT FIELD, 

4. TE INCIDENT FIELD, 



Computed Field dot Ideal Field 
Scaled Error component 
Computed Field dot Ideal Field 
Scaled Error component 



Indicate your selection entering "1","2","3", or "4" 
==> 2 

CALL PLOTTING ROUTINE 

Enter 0 for 3D-stick or 1 for 3D-fill option 

==> 0 

Default view angles are: phi=-45.,theta=70. 

(** SEE FIGURE 4.5 **) 
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The datafile used to generate fields is intfil.INT 
TEST DATA FOR USER EXAMPLE 
THIS HEADING FOLLOWS THE DATA THROUGHOUT 
Graphics caption 

TM INCIDENT FIELD, Scaled Error component 



0. 50 r- 




Inc. 



win. Z = +. 013; 
A = + . 000; Plane G 



wax. Z = 

= +.000; 



*.436; average Z =• +. Z27 
Uiew; phi=+80.0 theta=+70. 0 



Figure 4.5 3-D plot of scaled error component for example ogive radome 
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ENTER NEW ANGLES PHI & THETA (DEG) OR (9,9) WHEN DONE 



==> 9,9 

CORFLD OUTPUT GRAPHICS PRESENTATION 
CORE FIELD REPRESENTATION SELECTION MENU 

0. FINISHED WITH THIS ANGLE OF INCIDENCE 

Indicate your selection entering ’T","2",...,”10",or "0" 

==>0 

********* CORE ANALYSIS COMPLETED ************** 

Stop - Program terminated. 

When the user is finished with the angle of incidence selected, "0" is 
selected at the selection menu and the program execution is 
terminated. 
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V. PROGRAM VALIDATION AND RESULTS 



In order to validate the execution of EMRAD and CORFLD, the 
tangent ogive radome generated in Chapter IV was modified to have 
a relative permittivity equal to that of air, e r =1.0. The field 
coefficients were generated by EMRAD and the core fields were 
computed by CORFLD and compared to the unperturbed incident field 
as listed in Chapter III. 

Figures 5.1 and 5.2 show the magnitude and phase of the 
scaled dot product, as defined in Chapter III, between the computed 
and theoretical values for the TM and TE incident fields, respectively, 
for a "nose-on" incident angle. Each comparison shows a 94.9% 
magnitude matching with an almost constant 2.59 0 phase differential 
across the array. Figure 5.3 exhibits the scaled error component for 
the TM and TE incident fields, with an average magnitude of error 
component of 0.4%. The comparisons for this "nose-on" angle of 
incidence show that there is good agreement between the theoretical 
and calculated field components. It is observed that the TM and TE 
components displayed are 90 0 out of phase with each other. This is 
expected for "nose-on" incidence since the TM wave has its direction 
component in the -x direction, while the TE wave's direction 
component is in the y direction. 

Figure 5.4 shows the magnitude and phase of the scaled dot 
product between the computed and theoretical values for the TM 
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T lie datafile used to generate fields is airfil.INT 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 



TM INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 

1.50 



1. 00 




1.00 



min. Z = + .98Z; max. Z = +.996; average Z = +.949 
Inc. A = +.000; Plane G = +.000; Uiew: phi=+45.0 theta=+75. 0 



The datafile used to generate fields is airfil.INT 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 



TM INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 




min. Z = +1.73; max. Z = +2.95; average Z = +2.59 
Inc. A = +.000; Plane G = +.000; Uieu: phi=+80.0 theta= + 75. 0 



Figure 5.1 Dot product magnitude and phase for TM 
computed fields determined for radome composed of air 
compared with theoretical plane waves 
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The datafile used to generate fields is airfil.INT 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 



TE INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 




The datafile used to generate fields is airfil.INT 
Shaping function = TAN ogive 
Material of fornation is air, Er = 1 



TE INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 




Inc. A = +.000; Plane G = +.000; Uiewl phi=+80.0 theta= + 75.0 



Figure 5.2 Dot product magnitude and phase for TE 
computed fields determined for radome composed of air 
compared with theoretical plane waves 
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The datafile used to generate fields is airfil.INT 
Shaping function = TAN ogive 
Material of formation is air, Er = i 



TE INCIDENT FIELD, 
0.10 



Scaled Error component 




1. 00 



-0.05 c - 



min. Z = +.000; max. Z = +.023; average Z = +.004 
Inc. A = +.000; Plane G = +.000; Uieu! phi=+10. 0 theta= + 80. 0 



The datafile used to generate fields is airfil.INT 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 

TM INCIDENT FIELD, Scaled Error component 

0.10 r- 



0.05 — 




min. Z = +.000; max. Z = +.023; average Z = +.004 
Inc. A = +.000; Plane G = +.000; Uiew: phi=+75.0 theta=+80. 0 

Figure 5.3 Scaled error component for TM and TE 
computed fields determined for radome composed of air 
compared with theoretical plane waves 



76 



The datafile used to generate fields is airfil.INT 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 



TM INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 

3.00 



1. 00 




1.00 



win. Z = +.540; wax. Z = +1.66; average Z = +.874 
Inc. A = +60. 0; Plane G = +60.0; Uieu? phi=+45.0 theta=+75.0 

The datafile used to generate fields is airfil.INT 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 

TM INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 

78.00 



1.00 




1. 00 



min. Z = -27.9; max. Z = +60.1; average Z = -2.84 
Inc. A = +60. 0; Plane G = +60.0; View! phi = + 45. 0 theta= + 75. 0 

Figure 5.4 Dot product magnitude and phase for TM computed fields 
determined for radome composed of air compared with theoretical 
plane waves for incident angle of 60 degrees 
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incident field, for an incident angle 60° from the "nose-on” aspect. 
The magnitude of the scaled dot product between the computed and 
theoretical TM incident wave shows modal variation around the 
antenna plane. The magnitude of the computed value in the 
direction of the theoretical plane wave varies between 54.0% and 
166%, with an average magnitude of 87.4%. The phase component of 
the dot product varies across the plane from -27.9° to +60.1°, with 
an average phase of -2.84°. The maximum phase differential occurs 
at the edge of the antenna plane closest to the point of wave 
incidence, as the incident wave is defined as arriving in the -x - z 
axis. This variation seen in the magnitude and phase of the scaled 
dot product in the direction of the theoretical, undisturbed plane 
wave is confirmed in the scaled error component of Fig. 5.5 with an 
amplitude variation from 8.4% to 88.6%, averaging at 23.6% of the 
computed field evaluated as an error term perpendicular to the 
theoretical fields expected. 

Figures 5.6 and 5.7 exhibit the behavior of the calculated TE 
60° incident plane wave, which is similar, but less in magnitude 
compared to the calculated TM wave. Less modal variation is seen 
along the exhibited plane of the antenna. The scaled error 
component has a symmetry about the x-axis. The scaled error 
amplitude varies form 0.3% to 33.1%, with an average of 12.5%. 

The form of the error components as analyzed for this radome 
constructed of air identifies a phase error in the determination of the 
computed fields. The data computations of the EMRAD subroutines 
have been previously validated by Morgan. These validations 
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The datafile used to generate fields is airfil.INT 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 



TM INCIDENT FIELD, Scaled Error component 




min. Z = +.084# max. Z = +.886; average Z = +.236 
Inc. A = +60. 0; Plane G = +60.0; Uieu: phi=+45.0 theta=+75. 0 



Figure 5 5 Scaled error component for TM computed fields 
determined for radome composed of air compared with 
theoretical plane waves for incident angle of 60 degrees 
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The datafile used to generate fields is airfil.INI 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 



TE INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 




min. Z = *.860; max. Z = +1.34; average Z = +.981 
Inc. A = +60.0; Plane G = +60.0; Uiewl phi=+45.0 theta= + 75.0 



The datafile used to generate fields is airfil.INT 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 



TE INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 

27.00 — 




min. Z = -18.2; max. Z = +4.96; average Z = -.430 
Inc. A = +60.0; Plane G = +60.0; View*. phi=+45. 0 theta=+75.0 



Figure 5.6 Dot product magnitude and phase for TE computed fields 
determined for radome composed of air compared with theoretical 
plane waves for incident angle of 60 degrees 
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The dalafiie used to generate fields is airfii.INT 
Shaping function = TAN ogive 
Material of formation is air, Er = 1 

TE INCIDENT FIELD, Scaled Error component 

0.50 r- 




min. Z = +.003; max. Z = +.331; average Z = +.125 
Inc. A = +60.0; Plane G = +60.0; Uieu'. phi=+20. 0 theta= + 75.0 



Figure 5.7 Scaled error component for TM computed fields 
determined for radome composed of air compared with 
theoretical plane waves for incident angle of 60 degrees 
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considered the magnitude and phase of the scattered fields [Refs. 3 
and 4], and the magnitude of the core fields [Ref. 10] along the main 
coordinate axes. This work is the first to plot both the phase and 
magnitude variations across the entire two-dimensional plane and 
identifies the possible need for increased computational accuracy in 
the numerical routines of EMRAD. There is also the possibility of 
unresolved programming errors in CORFLD which should be further 
investigated. 

B. COMPUTATION RESULTS FOR VARIOUS RADOME MATERIALS 

EMRAD and CORFLD executions for the tangent ogive standard 
established in Chapter IV were run for varying material permittivity 
according to the material properties listed in Table 2.1. The 
comparisons of the computed fields to the unperturbed theoretical 
incident plane waves are included only for the "nose-on" incident 
angle which has good results for the ideal "air" construction. 
Fiber-reinforced glass, with e r =4.0; glass ceramic, with e r =5.5; and 
refractory oxide, with e r =8.0, are included. The form of field 
comparisons are similar in shape for each radome considered, only 
the magnitude of the plots varies. As expected from the discussion 
in Chapter II, the attenuation of the incident plane wave, received at 
the plane of the antenna in the core, increases with increasing 
relative permittivity. When the radome was constructed of air, the 
computed magnitude of each incident wave was 94.9%. When 

fiber-reinforced glass construction is considered, Figures 5.8, 5.9, and 
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5.10, the average magnitude of the computed wave is 60.1% and the 
average scaled error component is 22.7%. For the radome 
constructed of glass ceramic, Figures 5.11, 5.12, and 5.13, the average 
magnitude is 66.9% and the average scaled error component is 13.3%. 
For the refractory oxide radome, Figures 5.14, 5.15, and 5.16, the 
average magnitude of received field is 87.6% and the average scaled 
error component is 20.6%. 
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The datafile used to generate fields is P4FIL. 1NT 
Shaping function = TAN ogive 

Material of forwation is fiber-reinforced glass 



TM INCIDENT FIELD, Conputed Field dot Ideal Field, MAGN. 




win. Z = +.475# wax. Z = +.811; average Z = +.601 
Inc. A = +.000; Plane G = +.000; View: phi=+35.0 theta=+70. 0 



The datafile used to generate fields is p4fiI.lNT 
Shaping function = TAN ogive 

Material of forwation is fiber-reinforced glass 



TM INCIDENT FIELD, Conputed Field dot Ideal Field, PHASE 




or nin. Z = -71.1; wax. Z = -27.4; average Z = -51.8 



Inc. A = +.000; Plane G = +.000; View: phi=+10.0 theta= + 85. 0 

Figure 5 8 Dot product magnitude and phase for TM computed fields 
determined for radome composed of fiber-reinforced glass compared 
with theoretical plane waves for incident angle of 0 degrees 
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The datafile used to generate fields is p4fil.INI 
Shaping function = TAN ogive 

Material of formation is fiber-reinforced glass 




win. Z = +.475; wax. Z = +.811; average Z = +.601 
Inc, A = +.000; Plane G = +.000; View: phi=+45.0 theta= + 75.0 



The datafile used to generate fields is p4fil.INT 
Shaping function = TAN ogive 

Material of formation is fiber-reinforced glass 



TE INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 




Figure 5.9 Dot product magnitude and phase for TE computed fields 
determined for radome composed of fiber-reinforced glass compared 
with theoretical plane waves for incident angle of 0 degrees 
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The datafile used to generate fields is p4fil.INT 
Shaping function = TAN ogive 

Material of formation is fiber-reinforced glass 
TM INCIDENT FIELD, Scaled Error component 




win. Z = +.013; max. Z = *.436; average Z = *.227 
Inc. A = *.000; Plane G = +.000; View: phi=+95.0 theta=+70. 0 



The datafile used to generate fields is p4fil.INT 
Shaping function = TAN ogive 

Material of formation is fiber-reinforced glass 
TE INCIDENT FIELD, Scaled Error component 




min. Z = +.013; max. Z = +.436; average Z = +.227 
Inc. A = +.000; Plane G = +.000; View: phi=+10. 0 theta= + 70. 0 



Figure 5-10 Scaled error component for TM and TE computed fields 
determined for radome composed of fiber-reinforced glass compared 
with theoretical plane waves for incident angle of 0 degrees 



86 




Tlie datafile used to generate fields is c55fil.INT 
TEST DATA FOR THESIS 
Radone is tangential ogive, 
constructed of ceramic, Er = 5.5 

TM INCIDENT FIELD, Computed Field dot Ideal Field 



i 



MAGN. 




min. Z = +.668; max. Z = +.803; average Z = +.669 
Inc. A - +.000;Plane G = +.000; Uiew: phi=+45.0 tl>eta=+75.0 



Tlie datafile used to generate fields is c55fiI.INT 
TEST DATA FOR THESIS 
Radone is tangential ogive, 
constructed of ceramic, Er =5.5 

TM INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 




min. Z = -77.5; max. Z = -47.5; average Z = -61.9 
Inc. A = +.000;PIane G = +.000; UieuT phi=+45.0 theta=+75.0 

Figure 5.1 1 Dot product magnitude and phase for TM computed fields 
determined for radome composed of glass ceramic compared 
with theoretical plane waves for incident angle of 0 degrees 
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The datafile used to generate fields is c55fiI.INT 
TEST DATA FOR THESIS 
Radome is tangential ogive, 
constructed of ceramic, Er = 5.5 

TE INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 




min. Z = +.668; max. Z = +.803; average Z = +.669 
Inc. A = +.000;PIane G = +.000; Uieu: phi=+45.0 theta=+75.0 



The datafile used to generate fields is c55fil.INT 
TEST DATA FOR THESIS 
Radome is tangential ogive, 
constructed of ceramic, Er = 5.5 

TE INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 




min. Z = -77.5; max. Z = -47.5; average Z = -61.9 
Inc. A = *.000;Plane G = +.000; Uieu: phi=+45.0 theta=+75.0 



Figure 5.12 Dot product magnitude and phase for TE computed fields 
determined for radome composed of glass ceramic compared 
with theoretical plane waves for incident angle of 0 degrees 
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The datafile used to generate fields is c55fil.INI 
TEST DATA FOR THESIS 
Radone is tangential ogive, 
constructed of ceranic, Er = 5.S 

TM INCIDENT FIELD, Scaled Error conponent 




win. Z = ♦ . 005 ; nax. Z = +.312; average Z = ♦ . 133 
Inc. A = +.000;Plane G = ♦.000; Uiew: phi=*80.0 theta=*70.0 



The datafile used to generate fields is c55fil.INT 
TEST DATA FOR THESIS 
Radone is tangential ogive, 
constructed of ceranic, Er = 5.5 

TE INCIDENT FIELD, Scaled Error conponent 

0.50 r- 




nin. Z = +.005; nax. Z = +.312; average Z = +.133 
Inc. A = +.000; Plane G = +.000; Uieu'. phi=+10.0 theta= + 70. 0 

Figure 5.13 Scaled error component for TM and TE computed fields 
determined for radome composed of glass ceramic compared 
with theoretical plane waves for incident angle of 0 degrees 
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The datafile used to generate fields is 0X8FIL.1NT 
Shaping function = TAN ogive 

Material of fornation is refractory oxide# Er = 8 

TM INCIDENT FIELD, Computed Field dot Ideal Field, MAGN. 




min. 2 = +.634; max. 2 = +1.08; average 2 = +.876 
Inc. A = + .000 ;P lane G = +.000; Uiew: phi=+45.0 theta=+75.0 



The datafile used to generate fields is 0X8FIL.INT 
Shaping function = TAN ogive 

Material of formation is refractory oxide, Er = 8 



TM INCIDENT FIELD, Computed Field dot Ideal Field, PHASE 



1.00 



: - 0 . 




Q Of> 

U.0 ^ 



1.00 



-83.00 [ — 

min. 2 = -68.3; max. 2 = -39.4; average 2 = -52.4 
Inc. A = +.000;Plane G = +.000; Uiewl phi=+80.0 theta=+75.0 



Figure 514 Dot product magnitude and phase for TM computed fields 
determined for radome composed of refractory oxide compared 
with theoretical plane waves for incident angle of 0 degrees 
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The datafile used to generate fields is 0X8FIL.INT 
Shaping function = TAN ogive 

Material of formation is refractory oxide* Er = 8 



TE INCIDENT FIELD* Computed Field dot Ideal Field, MAGN . 

2.00 ft 



1.00 




0.50 



1.00 



min. Z = +.634; max. Z = +1.08; average Z = +.876 
Inc. A = +.000;Plane G = +.000; Uieu! phi=+55.0 theta=+75.0 



The datafile used to generate fields is 0X8FIL.INT 
Shaping function = TAN ogive 

Material of formation is refractory oxide* Er = 8 

TE INCIDENT FIELD, Computed Field dot Ideal Field* PHASE 




1.00 0.50 0.00 -0.50 7.00 

f mnzE 

i 



”0r| 

^01:00 



r*in. Z = —68.3; nax. Z = -39.4; average Z = -52.4 
Inc. A = * .080; Plane G = +.000; View: p}ii=+85.0 t!ieta=+80.0 



Figure 5.15 Dot product magnitude and phase for TE computed fields 
determined for radome composed of refractory oxide compared 
with theoretical plane waves for incident angle of 0 degrees 



91 




The datafile used to generate fields is 0X8FIL.INT 
Shaping function = TAN ogive 

Material of formation is refractory oxide# Er = 8 
TM INCIDENT FIELD, Scaled Error component 




min. Z = ♦ .006; max. Z = -*-.568; average Z = ♦.206 
Inc. A = ♦.880;Plane G = ♦.000; View: phi=+85.0 theta=+75.0 



The datafile used to generate fields is 0X8FIL.INT 
Shaping function = TAN ogive 

Material of formation is refractory oxide# Er = 8 



TE INCIDENT FIELD# 
1.50 



Scaled Error component 



1.00 




1.00 



min. Z = ♦.086; max. Z = ♦.560; average Z = +.206 
Inc. A = ♦ .000 ; Plane G = ♦.000; Uieu: phi=*25.0 theta=^75.0 



Figure 5.16 Scaled error component for TM and TE computed fields 
determined for radome composed of refractory oxide compared 
with theoretical plane waves for incident angle of 0 degrees 
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VI. CONCLUSIONS 



This thesis has presented a useful set of programs to 
numerically model a radome, designed of specified shape, 
construction, and material, and to evaluate the perturbation of the 
electromagnetic field as it passes through the modelled radome to 
the antenna within the radome structure. The radome design 
methodology was discussed in Chapter II, along with the radome 
shaping functions, wall construction, and choice of materials for 
radome construction. The program RADSHP is introduced and the 
shaping function, as well as the shaping factors for each function are 
explained. Chapter III developed the electromagnetic theory 
required to assemble the field components within the interior 
radome core, using the coefficients provided by the EMRAD program. 
The theory required for transforming the coordinates of the canted 
antenna plane to the core coordinates for an assumed two-axis 
gimbal was also discussed. Chapter IV explained the data flow for 
the radome model design process by using an example of the 
program executions for the design of a monolithic tangent ogive, with 
computational analysis by EMRAD and field component generation by 
CORFLD. The verification of the radome design and optimal placing of 
the core origin by the use of CURVE DIGITIZER software is also 
discussed. Chapter V discussed the computational errors discovered 
in the 3-D plotting of the field components. The scaled magnitude of 
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the computed perturbed plane waves for the "nose-on" incident 

angle agreed within 95% of the theoretical value and the phase of the 
computed field was an almost constant 2.59° out of phase with the 
theoretical incident plane wave. The errors increase when the angle 
of incidence is off the "nose-on" aspect, with modal amplitude 
variations from 54% to 166% of the theoretical plane wave. Further 
diagnostics are required to determine the source of this error. 

This is only the first stage in the development of a 

mathematical model for the interaction between the radome and the 
antenna system that the radome is designed to enclose. The model, 
in its present form, does not consider the presence of the antenna 

structure within the radome core. Further work is needed to 
determine the effect of the antenna structure and its interaction with 
the radome. When such a model is developed the predicted data can 
be compared to data collected by a physical antenna system. The 
culmination of this effort could be the analytic development of 
corrective lenses, as referred to in Chapter II, to correct or minimize 
the perturbation of the incident electromagnetic field caused by the 
radome structure. 
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Appendix A 



EMRAD HELD THEORY 

The program EMRAD uses the unimoment method mentioned in 
Chapter II, as developed by Morgan and Mei [Ref. 5] to decouple the 
interior and exterior spherical solutions. The expansion coefficients 
for the electromagnetic fields within the radome structure are 
generated as part of the scattering solution. CORFLD uses these 
coefficients for construction of the interior fields. The theory 
required by CORFLD to formulate the electromagnetic fields within 
the radome is discussed in Chapter II. 

The code developed by Morgan for EMRAD uses an optimized 
variational finite-element algorithm, in conjunction with a 

tri-regional unimoment method, to provide scattering solutions for 
each of multiple incident fields impinging on the specified structure. 
EMRAD represent a significant improvement in computational speed 
and versatility over computer codes which find the solution to the 
electromagnetic scattering and penetration problem by use of surface 
current integral equation formulations. This advantage results 
because of the highly sparse matrices embodied in the finite element 
solution vis-a-vis full matrices found using an integral equation 

formulation. An overview of the finite element scattering algorithm, 
along with the concept of coupled azimuthal potentials and the 

unimoment method are included in this appendix. Further 

information is available in references 3, 4, 5, 6, and 7. 
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A. COUPLED AZIMUTHAL POTENTIALS 

In the radome construction the material constitutive 
parameters, e(r) and p(r), are invariant to the azimuthal coordinate 
due to the axial symmetry of the radome. This axisymmetry allows 
implementation of vector field generation via coupled azimuthal 
potentials, as developed by Morgan, Chang and Mei [Ref. 7]. 

For ease of solution and representation, a normalized 
cylindrical coordinate system is adopted, as defined by (R , Z , <t>) = 
(k 0 p, k 0 z, <)>) where (p, z, ty) are standard circular cylindrical 
coordinates and k 0 = 2ttA 0 is the free-space wavenumber of the 

time-harmonic field. The equations in this form are related to the 
standard spherical coordinates (r, 0, <)>) as shown in Fig. A.l using the 
substitutions 



and 



R = k 0 r sin0 



(A.l) 



Z = k 0 r cos0 . (A. 2) 

The total electromagnetic fields are decomposed into azimuthal 
modes by an exponential Fourier series in the <]) -coordinate such that 



E(R,Z,<|>) = £ e m (R,Z) <>♦ 

m=— 

and 



(A. 3) 



ri 0 H(R,Z,<|>) = f h m (R,23 ei" 1 * , (A.4) 

m=-«> 
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X 



Figure A.l Rectangular (x, y, z), spherical (r,0, <t>), and 
normalized cylindrical (R, Z,8) coordinate systems 
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where T| 0 = 120 7t Q is used to normalize the magnetic fields, so that 
h = tj 0 H. 

Using the axial symmetry of the radome construction, the 
modal fields decouple when the field expansions are substituted into 
Maxwell's equations, yielding [Ref. 7], after judicious cross 
substitutions the following modal field generating equations in terms 



of two scalar potentials, T]( R,Z,m ) and T 2 ( R,Z,m ): 

4> x e m = j 5n ( m<J>xVrj - R^ vr 2 ) (A. 5) 

<f • e m - £ (A.6) 

<l>xhm = j f m ( m<>xVr 2 + Re r VIj ) (A. 7) 

♦ • h m = S. (A.8) 

where the two dimensional gradient operator is defined as 

V = [ R (9/3r ) + Z (9/9Z) ] (A. 9) 

and f m is a multiplicative function given by 

UR.Z) = [ n r (R,2) e,(R,Z) R 2 - m 2 J 4 . (A.l 0) 

Note that the potentials, T] and r 2 are proportional to the azimuthal 



components of the electric and magnetic fields, e m and h m . These 

potential functions are continuous everywhere, including at dielectric 
and magnetic interfaces. This property of uniform continuity in the 
field is used to obtain a numerical solution of the fields. 
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The total fields may be described in term of the uniformly 
continuous scalar potentials Ti( R,Z,m ) and T 2 ( R,Z,m ) such that 

e m (R,Z) = j f m (mVli - Rn, ♦ x Vr 2 ) + * (A.l 1) 

and 

h m (R,Z) = j f m (m Vr 2 + Re r <|> x VT, ) + + ^ . (A.l 2) 

This method of expanding the electromagnetic fields into modal 
components is named the Coupled Azimuthal Potential (CAP) method 
due to the fact that the scalar potentials Ti(R, Z, m) and T 2 (R, Z, m) 
are coupled in these field equations and are not separable. 

Due to the symmetry in the azimuthal coordinate, <)>, in radome 
design, an arbitrary meridional cross section of the radome 
construction may be represented in an (R, Z)-plane such that there 
exists a two-dimensional solution domain 5, with a boundary curve 
forming the surface of revolution 55. It has been shown by Morgan, 
Chang and Mei [Ref. 7] that in this solution domain the potentials, V x 
and T 2 , satisfy the following formally self-adjoint pair of partial 
differential equations (PDE's), 

v-(f m Re r vr,)+-|-r 1 + mV(M*vr 2 ) = o (A.l 3) 

and 

-mV-t^xVIi) + V-(k R| 1 , vr 2 )+-|-r 2 = 0 , (A.l 4) 

with Dirichlet boundary conditions (BC) for Tj and r 2 specified on 55 . 
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In addition to this PDE formulation there is an easily managed 
variational principle [Ref. 7] established using either an 
Euler-Lagrange approach or a stationary-theorem approach. The 
variational principle is based upon the existence of a stationary 
functional of the potentials, T j and T 2 , and their gradients such that 

sf = || L(R,zr 1 ,r 2 ,vr„vr 2 ) dRdz , (A.i5) 

S 

where the Lagrangian operator, L, is unique within an arbitrary 

constant multiplier and arbitrary independent additive function. The 
Lagrangian operator has the quadratic form 

L = u V 15 . (Re r Vrj + m<|.xvr 2 ) 

+ ^ vr 2 - (Rn r vr 2 - m$x vrj ) 

- ^ (e r if + Hr r 2 2 ) . (A.16) 

With this analytical approach to the solution of the electromagnetic 
field's behavior through the axially symmetric radome structure, the 
determination of the field on a point-to-point basis is accomplished 
through a numerical finite-element algorithm. 

B. FINITE ELEMENT ALGORITHM 

The numerical solution of the electromagnetic field problem 

involving scattering by, and penetration of, the radome is 

accomplished by the development of a finite-element algorithm 

using the mesh geometry in Fig A. 2. This configuration in the 

meridional (R, Z)-plane uses linear triangular elements which are 

1 00 




Unbounded 
Exterior 
REGION III 



Figure A. 2 Finite Element Mesh Geometry 
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arranged by the EMRAD program to conform to the surfaces of the 
axisymmetric radome. The selection of surface conformation 
simplifies numerical analysis of the normal derivative discontinuities 
at the surface boundaries. The two-dimensional mesh configuration 
that is developed is named the semiannular conformal (SAC) mesh. 

Figure A. 2 graphically represents the separate spatial regions 
which are of interest. The homogeneous inner spherical core of the 
radome occupies Region I. The inhomogeneous axisymmetric radome 
wall construction cross-section is in Region II. The unbounded 

exterior region exists exterior to the outer sphere in the 
homogeneous free-space of Region III. Note that the surfaces of the 
inner sphere, Region I, and the exterior Region III, bound the area of 
interest for numerical analysis in Region II. 

The electromagnetic fields within the homogeneous inner core 
of the radome (Region I) and the unbounded exterior (Region III) are 
defined in terms of the appropriate spherical harmonic modal 
expansions of the radial TE and TM field components. At the 
surfaces, r = a and r = b, the modal elements in these fields are 
applied as sets of Dirichlet boundary conditions. Then, the numerical 
solution for the fields in the finite-element domain of Region II is 
generated for each boundary condition. Finally, the unknown modal 
expansion coefficients for the core and exterior scattered fields are 
obtained by enforcing continuity of the electromagnetic fields at the 
interfaces along r = rj and r = r 2 . This approach to numerical solution 

of the fields is the basis of the unimoment method. In the 
implementation of the variational principle using the finite-element 

1 02 



method, the azimuthal modal field components are selected as the 
fundamental unknowns. 

In Region II of Fig. A. 2, the SAC mesh is composed of numerous 
interior nodes. Each of these interior nodes is surrounded by six 
triangular regions known as "elements". The azimuthal field 
components within the mesh are represented as basis function 
expansions, 

N 

♦ • e m (R,Z) = £ e m (n) u„(RZ) (A.17) 

n=l 

and 

♦ • h m (R,Z) = £ h m (n) u„(R2) . (A.l 8) 

n=l 

Each of the respective sets of e m (n) and h m (n) represent one of the N 
complex nodal values found in the azimuthal fields. The finite 
element basis functions u n (R, Z) are area coordinate linear pyramid 

functions which each assume a value of unity at their associated 
node, n, and zero at each of the neighboring nodes. There is an 
assumption of linear variation within each of the triangular elements 
in the surrounding group such that 

u n (R,Z) = a n> ,R + b M Z+ c M . (A.l 9) 

The real constants (a, b, c) corresponding to each node, n, in element 

£ is found in terms of the nodal (R, Z)-coordinates of the element. 

The variational finite-element technique is initiated by 

substituting the linear shape-function expansions in Equations (A.17) 

and (A. 18) into the stationary Euler-Lagrange functional, SF, as 
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described in Equation (A. 15). The stationary point of SF is then 
computed by imposing the null derivative condition with respect to 
each of the unknown coordinates representing the interior nodal 
fields, while constraining the "known" boundary nodal fields. The 
resulting linear system describing the nodal fields is 

N 

S [f m Re r V(Ru n ) V(Ru i )-R£ r u„u i ]dRdZ 

i=1 S n 

+h m (i) jj [m^V(Ru n ) *xV(R Ui )]dRdZ)=0 (A.20) 

s n 

and 



I{h m (i)JJ [(„ R 1 i I V(Ru„)-V(Ru i )-Rn I u n u i ]dRdZ 
w Sn 

-e m (i)JJ [mt,V(Ru n ). + xV(Ru|)]dRdZ}=0 (A.2 1) 

Sn 

for each value of the internal nodes, n, where i spans all nodes, 
including boundary nodes. As these above equations are enforced 
for each internal node, they involve integrations over only the 
elements connected to this node, S n , and therefore relate the value of 

e m (n) and h m (n) only to adjacent nodal values, e m (i) and h m (i). Then, 

in the SAC mesh implementation, in the case of each internal node 
and the six adjacent nodes, the above equation will relate 14 
adjacent node azimuthal field values, some of which may be specified 
at the boundary. 
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The computation for each nodal value, as indicated by 
Equations (A.20) and (A.21), results in a sparse global system matrix 
which contains generally complex elements. At each node the mean 
values of the spatially variable e r (R, Z) and p r (R, Z) are used as the 
material constants in that triangular element. Applying the Dirichlet 
boundary conditions to the SAC mesh shown in Fig. A. 2, Morgan and 
Mei [Ref. 5] found the numerical finite-element solution for the 
electromagnetic fields in Region II using a two-sweep block-by-block 
matrix inversion algorithm related to the Riccati transform for 
difference equations. At each constant colatitude, 0 = 0 k , a vector of 
internal nodal azimuthal field values along the ray from rj to r 2 is 
defined by p k • The resultant global system matrix defined by 
Equations (A.20) and (A.21) has a diagonal tri-block submatrix 
structure, which can be represented locally by relationships between 
adjacent radial node vectors in the mesh, 

Ak Pk-i + Bk Pk + Q Pk+i = a k » (A. 2 2) 

such that A k , B k , and C k are sparse submatrices. These form the 
tri-block diagonal structure of the global system matrix while a k , 
the local driving function, is generated from Equations (A.20) and 
(A.21) using the known Dirichlet BC's at r=a and r=b. 

The system solution in the two-dimensional meridional (R, Z) 

plane is performed in a step-wise manner beginning at the positive 

Z-axis (0 k = 0j =0) and proceeding in equal marching increments of 

A0. At each 0 k , the SAC mesh is constructed using a minimum 

distortion scheme to conform to the approximate piecewise smooth 
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surfaces. At each 0 k increment, as the mesh is constructed in this 
sequential circumferential manner, the three block matrices are 
generated which relate adjacent radial node vectors. 

During the incremental marching of 0 k , the block-by-block 
Riccati matrix inversion is implemented. The Riccati matrix, R k , and 
the associated vector, S k , are related to the vector, P k , by the 
transformation 

Pk-, = Rk Pk + S k . (A. 2 3) 

Substituting this Riccati relation into Equation (A. 22), solving for P k 
and returning to the form of the Riccati transform yields the 
recurrence relations 

R k+ i = -( Bjj + A k »R k ) 1# Q (A. 24) 

and 

$k+i = ( + A k *R k ) 1# ( Pk ~ A k *S k ) • (A. 2 5) 

Using these relations, the Riccati transform sweeps forward 
from initial conditions, using the equations of step-wise evolution in 
Equations (A. 24) and (A.25) to sequentially generate and store the R k 

matrices and the S k vectors. This is called the "first sweep". The 
solution vectors P k are then obtained by using the transform in 
Equation (A. 23) to perform a "backsweep" with calculations using the 
stored Riccati arrays. This sequentially localized technique for 
handling the Dirichlet boundary value problem blends easily with 
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} 



the unimoment method which is used to solve the electromagnetic 
field in Regions I and III of Fig. A.2. 

C. UNIMOMENT METHOD 

The unimoment method is a computational procedure which 
provides for the numerical solution of the incident, scattered and 
spherical core electromagnetic field expansions when dealing with 
the constructed SAC finite-element mesh of Fig. A. 2 The numerical 
solution is generated by representing the core (r<rj) and the exterior 

scattered (r>r 2 ) fields in terms of radial TE and TM spherical 
harmonic expansions, each with unknown core and scattered modal 
field coefficients. In order to solve for these unknown coefficients, 
the core expansion of Region I is applied as a boundary condition 
along r = a, while the boundary condition at r = b is determined as 
the addition of the incident field and the scattered field expansion in 
Region III. 

Superposition is then applied to the boundary conditions, 

composed of the core field expansion and the exterior field expansion 

in order to find a complete solution. The complete field solution is 

found from the addition of the partial field solutions generated from 

two sets of applied BC's. The first set of boundary conditions applied 

is composed of the modal spherical harmonic core fields along r = a 

with zero fields along r = b. The second set of boundary conditions 

contains the spherical harmonic exterior scattered and incident fields 

at r = b with null fields at r = a. These multiple boundary conditions 

generate multiple solutions at the nodes along both boundaries, 
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r = Tj and r = r 2 . The unknown core and scattered modal field 

coefficients are obtained by setting the weighted superposition of 
interior solutions equal to the original field expansions. This equality 
is enforced along the boundaries at r = r i and r = r 2 , utilizing least 
squares. 

The success of the unimoment method is contingent upon the 
capability to analytically represent the electromagnetic fields 
everywhere outside of the finite element mesh region. These 
representations, as developed for Regions I and III, are given as 

infinite series of spherical harmonics which must be truncated to 
facilitate numerical calculations. The plane-wave fields that are 

incident on the radome are decomposed in spherical coordinates into 
azimuthal modes as 



M 



El<r.<M>) = I <4,^,6) 



m=0 



and 



M 



HoHKrAO) =1 h^(r,0) 



m=0 



jsin(m<t)) 
cos(m<{>) 

cos(m<j)) 
jsin (m<J))_ 



(A. 2 6) 



(A. 27) 



where M is the mathematical truncation index of the series 
representation. The bracketed sinusoids represent the polarization 
convention used for TE and TM modes, respectively. 

During computations, the numerical solution of the unknown 



core and scattered modal field coefficients is run concurrently with 

the sequential Riccati transform finite-element solution in Region II. 
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On the forward sweep, boundary conditions for both the incident and 
the modal fields are generated only at the local node row being 
considered in that particular sweep step. On the backsweep the 
nodal values for each applied boundary condition are obtained along 
the r = Tj and r = r 2 contours. The corresponding incident field and 

spherical harmonic field at each node are generated, and then the 
nodal residual fields are formed. During the backsweep, inner 
products integrations are accumulated and their contributions are 
added to a field moment matrix and the incident-field driving vector, 
which corresponds to each incident field. The field expansion 
coefficients for each incident field are obtained by multiplying the 
inverse of the field moment matrix with that field's driving vector. 
With the field expansion coefficients calculated and collected, the 
electromagnetic field within the radome core can be determined on a 
point-to-point basis as developed in Chapter II . 
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RADSHP PROGRAM CODE 



PROGRAM RADSHP 
C 

£ ************************************************************** 

C 

C GENERATING THE MATERIAL FILE AND DATA STRUCTURE FILE 
C FOR A LAYERED DIELECTRIC RADOME 
C MATERIAL FILE PREPARED FOR EMRAD 
C DATAFILE PREPARED FOR CURVE DIGITIZER 
C USE EMCADIN TO GENERATE STRUCTURE FILE FOR EMRAD 

C Written by R.J. Vince Nov-Dec 89 
C 

£ ************************************************************** 

C VARIABLE DECLARATIONS 

q *************************************************************** 
REAL X(5,1001), Z(5,1001), PI, DELTHETA, DEGTORAD 
REALTHETA,C(5),R(5),H(5),DELZ(5),SHFTEACTORE(5,7) 

INTEGER IBIG, I, L, SEPRATN, NSHAPE, NLMAX, RESOLV, RESPTS 
COMPLEX ER(5),UR(5) 

CHARACTER* 10 SHAPE(6) 

CHARACTER * 64 HDR1.HDR2 
CHARACTER* 8 YNAME 
CHARACTER* 12 MATFIL.DATFIL 
CHARACTER*! YN 



SLARGE 

C 

Q ******** ************************************************ ******* 

C INITIAL VALUES 

p *************************************************************** 
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PI=3.14 15927 
DEGTORAD = PI/180. 

SEPRATN = 999990 
RESOL V= 1000 
RESPTS=RESOLV+ 1 
NLMAX = 5 

C ITITIALIZE SHAPE FUNCTION FACTOR ARRAY 
C 

DO 67 1=1 .NLMAX 
DO 65 J=l,6 
F(U)=0.0 
65 CONTINUE 
67 CONTINUE 
C 

C 

Q ****************************************************************** 

C INPUT NECESSARY DATA 

£ ****************************************************************** 

C 

C PROMPT FOR OUTPUT FILENAMES 
C 

WRITE(MIO) 

110 FORMAT(//9X,'WELCOME TO RADSHP’, 

1 ///9X.THIS PROGRAM GENERATES RADOME MATERIAL AND SHAPE DATA', 

2 //7X/FILES FOR USE WITH EMRAD, BASED ON USER INPUT.’, 

3 ////9X.TWO FILES ARE REQUIRED TO HOLD PROGRAM OUPUT DATA.', 

4 ///9X; PLEASE ENTER THE NAME OF THE FILE TO BE USED TO 

5 /7X/HOLD THE PROGRAM MATERIAL PARAMETER DATA. THE EXTENSION’, 

6 /7X/.MAT WILL BE ADDED AUTOMATICALLY, I. E. YOURNAME.MAT’,/) 
READ(*,100) YNAME 

NCHAR = INDEX(YNAME,' •) -1 
MATFIL = YNAME(1:NCHAR) I I '.MAT 

OPEN( 1 TlLE=MATFIL,STATUS=’UNKNOWN’,ACCESS=’SEQUENTIAL', 

1 1 1 



1 FORM='FORMATTED) 

WRITER, 111) 

1 1 1 F0RMAT(///9X,’PLEASE ENTER THE NAME OF THE FILE TO BE USED 

1 /7X,TO HOLD THE PROGRAM OUTPUT DATA. THE EXTENSION .DAT WILL', 

2 /7X,WILL BE ADDED AUTOMATICALLY, I. E. YOURNAME.DAT /7X) 
READ(MOO) YNAME 

NCHAR = INDEX(YNAME,' 0 -I 
DATFTL = YNAME(1 :NCHAR) // '.DAT 

OPEN(2,FTLE=DATFTL,STATTJS='UNKNOWN'ACCESS='SEQUENTIAL', 

1 FORM='FORMATTED') 



C 

C PROMPT FOR # OF LAYERS WITHIN THE RADOME 
C 

5 WRITE(*,114)NLMAX 

1 14 FORMAT(//7X,'ENTER THE NUMBER OF LAYERS (.LE. ',I2,’)7/7X) 
READ(*,*) NL 

IF ( (NL .LT. 1) .OR. (NL .GT. NLMAX) ) THEN 
WRITE(* , 1 1 5)NLMAX+ 1 

1 15 FORMAT(//7X,'NUMBER OF LAYERS IS OUT OF RANGE, ((k#<U2,')7/7X) 
GOTO 5 

END IF 



C LOOP TO INPUT LAYER INFORMATION 
C 

WRITE(*,7) 

7 FORMAT(/8X,'Each layer of the radome will have material consts,', 

1 /10X,'an outer height of the radome layer, H,', 

2 /lOX.’a outer base radius, R, and', 

3 /10X,'a shape function, & shape factor parameters (if appl.),', 

4 //lOX.'The First "layer" will be the inner core of the radome, ’, 

5 /lOX.'normally filled with air.') 
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39 



DO 12 L=1,NL 
CONTINUE 
WRITE)* ,9) L 
9 FORMAT)/,’ 1 Layer No: ’43) 

IF (L.EQ.l) THEN 

WRITE)*,*) ' This is the inner core of the radome' 
ELSEIF (L.EQ.NL) THEN 

WRITE)*,*) ' This is the outer layer of the radome' 
END IF 



C Prompt for Height and Radius of this layer 
C 

WRITE)*,*) 

WRITE(*,43) 

43 FORMAT(/8X, Please enter the radome layer height, H,' 

1 /8X,' and the outer base radiusjt, for this layer.') 

RE AD(*,*) H(L)P(L) 

C Prompt for shaping function, to determine equations to be used 
C to generate this layer of the radome 
C 

WRITE)*,*) 

WRITE)*, 10) 

10 FORMAT(/8X,’The following shaping functions are available, ' 

1 /8X,' 1. General Ogive' 

2 /8X,' 2. Tangent Ogive' 

3 /8X,’ 3. von Karman' 

4 /8X,' 4. Power Series' 

5 /8X,' 5. Parabola, curved' 

6 /8X,' 6. Parabola, capped by pointed apex' 

7 //8X,' Indicate desired shaping function by entering a number' 

8 /10X,' 1,2, ...or 6’) 



1 1 3 



READ(V) NSHAPE 



IF (NSHAPE.EQ.1) THEN 
C Ogive 

SHAPE(L) = 'Gen. Ogive' 

CALL OGIVE(L^^,H,R ( F,RESOLV > RESPTS^ILMAX) 

ELSEIF (NSHAPE.EQ.2) THEN 
C Tangent Ogive 

SHAPE(L) = TAN Ogive ' 

CALL TANOGIVE(L,X,Z,H > R ,F,RESOLV .RESPTS .NLMAX) 

ELSEIF (NSHAPE.EQ.3) THEN 
C von Karman 

SHAPE(L) = 'von Karman' 

CALL VONKRMN(L,X,Z,H,R,RESOLV,RESPTS,NLMAX) 

ELSEIF (NSHAPE.EQ.4) THEN 
C Power Series 

SHAPE(L) = "Pwr Series’ 

CALL POWER(L,X,Z,H^F^ESOLVRESPTS > NLMAX) 

ELSEIF (NSHAPE.EQ.5) THEN 
C Parabola (curved) 

SHAPE(L) = 'Para Curve' 

CALL PARACRV(L,X,Z,H,R,F,RESOLV,RESPTS .NLMAX) 

ELSEIF (NSHAPE.EQ.6) THEN 
C Parabola with apex 

SHAPE(L) = 'Para APEX' 

CALL APEX(L,X,Z,H,RF,RESOLV,RESPTS > NLMAX) 
ELSE 

C Selection out of limits 

WRITE (*,*) 'Selection is not available, try again.' 

WRITE (*,*) 
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GOTO 39 
ENDIF 



C Provide listing of inner layer parameters 
C 

WRITE(*,554) 

DO 11 1=1 JL 

WRITE (*,555) I,SHAPE(I),H(I)d?(I),(F(IT)J=l,7) 

11 CONTINUE 

WRITE(* *) 

12 CONTINUE 



C 

C SHIFT RADOME SO THAT ORIGIN IS NEAR THE CENTER OF INNER LAYER 
C NOTE: FACTOR is adjusted for optimal placement of origin for 
C mesh construction by EMRAD 
C 

FACTOR=0.4 



C Max inner height = Z(NLRESPTS) 

C 

SHFT=Z(1 1 RESPTS)*(F ACTOR) 

DO 88 L=1,NL 
DO 77 I=1,RESPTS 
Z(L,I)=Z(LJ)-SHFT 
77 CONTINUE 

88 CONTINUE 



C 

Q ****************************************************************** 

C OUTPUT 

£ ****************************************************************** 

1 1 5 



DO 330 L=1,NL 

C OUTPUT LAYER X.Z POINTS FROM TIP TO X-Z PLANE 
DO 220 I=RESPTS,1,-1 
WRITE(2 ,*) X(L,I),Z(L,I) 

220 CONTINUE 

C CLOSE LAYER BELOW SHIFTED X-Y PLANE 
C AT Z=[-1*[SUM OF (LAYER THICKNESS)] -SHFT] 

C 

C 

C DETERMINE LAYER THICKNESS, INNER LAYER HAS NO SHIFT FOR THICKNESS 

C OTHER LAYERS HAVE BOTTOM THICKNESS DETERMINED BY RADIAL THICKNESS 
C NOTE - THE SHIFT PERFORMED FOR THICKNESS IS ADDITIVE FOR LAYERS 
C 

IF (L.EQ.l) THEN 
THICK = 0 
ELSE 

THICK = THICK + R(L)-R(L-1) 

ENDIF 

C 

BOTTM = (-1.0 * THICK) - SHFT 
WRITE(2,*) R(L), BOTTM 
WRITE(2,*) 0.0, BOTTM 
WRITE(2,*) SEPRATN.SEPRATN 
330 CONTINUE 

C SIGNAL TERMINATION OF DATA SET 
WRITE(2,100) END, END' 



C Put summary of data to screen 
C 

WRITER ,5 54) 

DO 440 I=I,NL 

WRITE (*,555) I,SHAPE(I),H(I)J?(I),(F(I,J)J=1,7) 
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440 CONTINUE 



C 

C PROMPT FOR HEADER INFORMATION 
C 

WRITER, 112) 

1 12 FORMAT(//7X,THE HEADERS ALLOW THE USER TO IDENTIFY THIS7/7X, 
l'SET OF DATA FROM ALL OTHER SETS.’, 

2//7X/PLEASE ENTER HEADER #1 (64 CHARACTERS MAX)7/7X) 
READ(MOO) HDR1 
WRITER, 113) 

1 13 FORMAT(//7X,TLEASE ENTER HEADER #2 (64 CHARACTERS MAX)',//7X) 
READ(*,100) HDR2 

C 

C WRITE HEADER INFORMATION INTO MATERIAL OUTPUT DATA FILE 
C 

WRITE( 1,105) HDR1 
WRITE( 1,105) HDR2 

C 

C WRITE LAYER INFORMATION INTO MATERIAL OUTPUT DATA FILE 
C 

WRITE(1,106) NL 

C Prompt for material Er for each layer of construction 
C 

DO 812 L=1J4L 
WRITER, 9) L 

IF (L.EQ.l) THEN 

WRITE(*,*) ' This is the inner core of the radome’ 

ELSEIF (L.EQ.NL) THEN 

WRITE(*,*) ’ This is the outer layer of the radome' 

END IF 
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C Prompt for material constants 
C 

WRITE(*,*) 'Enter Er for this Layer: ' 
READ(V) ERR 

C Ensure that Imag Er is < 0 
C 

ERI = -0.000001 



ER(L)=CMPLX(ERR,ERI) 
UR(L)=CMPLX( 1 .0,- 1 £-6) 
WRITE(1,102) ER(L),UR(L) 
812 CONTINUE 



C PLACE COPY OF CONSTRUCTION AT BOTTOM OF MATERIAL FILE 
C 

WRITE(1,554) 

DO 660 1=1 ,NL 

WRITE (1,555) I,SHAPE(I),H(I)£(I),(F(U) r J=l,7) 

660 CONTINUE 



C — FINAL OUTPUT TO SCREEN 

WRITER ,999) MATFIL,DATFIL,DATFIL 
999 FORMAT(//8X, 

1' Computations and output are completed, the datafile ,',A, 
2/,8X,' holds the material specifications for input to EMCAD.' 
3//,8X,’ The datafile, ',A,', holds the radome shape data for' 

4/,8X,' viewing in the CAD package, CURVE DIGITIZER. ' 
5//,8X,' The structure file required for input to EMCAD can be' 
6/,8X,' generated by the executing the program EMCADIN,' 
7/,8X,' using, ’,A,' , as input to EMCADIN, ’ 

8//,8X,' The output file from EMCADIN will be in proper format' 
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9/,8X,' for input to EMCAD for the structure file.') 

WRITER, *) 

C 

C ****************************************************************** 

100 FORMAT(A) 

101 FORMAT(I5) 

102 FORMAT(4(E14.6)) 

103 FORMAT(I5,5(E14.6)) 

104 FORMAT(2(I5)) 

105 FORMAT(' » 

106 FORMAT(I5) 

554 FORMAT (///////////////// 

1 //1X,LAYER’,2X,’ SHAPE \2X,’HEIGHT,2X,' RADIUS', 

1 2X,' FI ',2X,’ F2 \2X,’ F3 \2X,’ F4 \2X/ F5 

2 2X,' F6 ',2X,' F7 VI X,'- — ',2X,' \2X,' ’, 

3 2x; — 'ax: — '.2X; — , ,2x;-- , ^x; — \ 

4 2 x; — \ 2 x: — \ 2 x ; — •) 

555 FORMAT (1X,I5,2X,A10,2X,F6.2,2X,F6.2,7(F6.2)) 

9999 CONTINUE 
C 

C END PROGRAM RADOME 
C 

CLOSE(l) 

CLOSE(2) 

STOP 

END 

*********************************** 

C OGIVE SHAPING FUNCTION SUBROUTINE 

SUBROUTINE OGIVE(L,X,Z,H,R,F,RESOLV,RESPTS,NLMAX) 

INTEGER IBIG, I, L, SEPRATN, NSHAPE, NLMAX, RESPTS, RESOLV 
REAL X(NLMAX,RESPTS), Z(NLMAXRESPTS), PI, DELTHETA, DEGTORAD 
REAL THETA,C(5),R(5),H(5),DELZ(5),SHFT .FACTOR 
REAL ALPHA 1(5), ALPHA(5), B, M, RL, GAMMA, F(NLMAX,7) 

COMPLEX ER,UR 



1 1 9 



CHARACTER*! YN 



C SET CONSTANTS 
C 

PI=3.1415927 
DEGTORAD = PI/180. 

C 

C DEFINE START AND END POINT OF EACH RADOME LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RESPTS)=H(L) 

X(L,RESPTS)=0.0 



C DEFINE RADOME LAYER BY OGIVE EQUATIONS .... 

DELZ(L)=H(L)/RES OLV 

C Prompt for Shaping factors for this ogive layer 
C 

WRITE(* ,*) 

WRITE(*,43) 

43 FORMAT(8X,'The ogive shaping function requires three shaping', 
1/8X, 'factors for determining the radome construction.' 
l//8X,'The first factor, FI, determine the arc of curvature', 
l/8X,’of the radome surface. Please enter FI:') 

READ(*,*)F(L,1) 

WRITE(*,*) 

WRITE(*,45) 

45 FORMAT(8X,' Two other shaping factors are used to determine the', 
1/8X, 'relative position of the center of the ogive shaping', 
l/8X,’arc, relative to the radome apex.', 

1//8X,’ The first of these factors, F2, determines the perp. ', 

1/8X, 'distance from the z-axis to the center of the arc.', 
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1//8X,' The second positioning factor, F3, determines the', 
l/8X,'the distance from the center of the arc to the apex', 
l/8X,'of the radome, along the z-axis.’, 

1//8X/Please enter the position factors for the center of the', 

1/8X,' curvature arc, F2 and F3. 1 ) 

READ{*,*)F(L,2) ( F(L,3) 

DO 33 I=2,RESPTS-1 

Z(L,I)=Z(L,I - 1 )+DELZ(L) 

C EQUATION OF RADOME IS x=SQRT( Fl**2-(F2-(H-z))**2 ) - F3 
C 

X(L4) = SQRT( F(L,1)**2 -( F(L,2) - (H(L)-Z(LJ)) )**2 )-F(L,3) 

33 CONTINUE 



C END OF RADOME LAYER .... 

RETURN 

END 

*************4^*******£f^D Qp QQJYP 
*********************************** 

C TANGENT OGIVE SHAPING FUNCTION SUBROUTINE 

SUBROUTINE TANOGIVE(L,X,Z,H,R,F,RESOLV,RESPTS,NLMAX) 

INTEGER IBIG, I, L, SEPRATN, NSHAPE, NLMAX, RESPTS, RESOLV 
REAL X(NLMAX, RESPTS), Z(NLMAX,RESPTS), PI, DELTHETA, DEGTORAD 
REAL THETA, C(5),R(5),H(5),DELZ(5),SHFT .FACTOR 
REAL ALPHA 1(5), ALPHA(5), B, M, RL, GAMMA, F(NLMAX,7) 

COMPLEX ER,UR 
CHARACTER* 1 YN 

C SET CONSTANTS 
C 

PI=3. 1415927 
DEGTORAD = PI/180. 
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c 

C DEFINE START AND END POINT OF EACH RADOME LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RESPTS)=H(L) 

X(L,RESPTS)=0.0 

C DEFINE RADOME LAYER BY TANGENT OGIVE EQUATIONS .... 

DELZ(L)=H(L)/RESOLV 

C Determine Shaping factors for this ogive layer 
C 

F(L,4) = R(L)/2.0 + 2.0*( H(L)**2.0 )/R(L) 

F(L,5) = F(L,4) - R(L)/4.0 

WRITE(*,*) 

WRITE(*,43)F(L,4),F(L,5) 

43 FORMAT(8X,'The tangent ogive shaping function determines two’, 

1/8X, 'factors for the radome construction. These are calculated’ 
l/8X,'based on the height and radius of the radome layer.’, 

1/8X, These factors are F4 = ’,F6.2,'and F5 = \F6.2) 

DO 33 I=2,RESPTS-1 

Z(L,I)=Z(L,I-1)+DELZ(L) 

C EQUATION OF RADOME IS x=SQRT( F4**2 - z**2 ) - F5 
C 

X(L,I) = 4.0* (SQRT( F(L,4)**2- Z(L,I)**2 ) - F(L,5)) 

33 CONTINUE 



C END OF RADOME LAYER .... 
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RETURN 

END 

**********************£Pjp) QP •p^pj OGIVE 
*********************************** 

C von Karman SHAPING FUNCTION SUBROUTINE 

SUBROUTINE VONKRMN(L,X,Z,H,R,RESOLV .RESPTS .NLMAX) 

INTEGER IBIG, I, L, SEPRATN, NSHAFii, NLMAX, RESPTS, RESOLV 

REAL X (NLMAX .RESPTS), Z(NLMAX JRESPTS), PI, DELTHETA, DEGTORAD 

REAL THETA,C(5),R(5),H(5)DELZ(5),SHFT,FACTOR,SQP,ZF 

REAL ALPHA 1(5), ALPHA(5), B, M, RL, GAMMA 

COMPLEX ER.UR 

CHARACTER* 1 YN 

C SET CONSTANTS 
C 

PI=3.14 15927 
SQP=SQRT(PI) 

DEGTORAD = PI/180. 

C 

C DEFINE START AND END POINT OF EACH RADOME LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RESPTS)=H(L) 

X(L,RESPTS)=0.0 

C DEFINE RADOME LAYER BY VON KARMAN EQUATIONS .... 

DELZ(L)=H(L)/RESOLV 

DO 33 I=2,RESPTS-1 

Z(L,I)=Z(L ,1- 1)+DELZ(L) 

C EQUATION OF RADOME IS 
C x=( R/(SQRT(PI)) ) * SQRT( FZ - SIN(2*FZ)/2 ) 
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C for FZ = ACOS( 1 - 2 * (H-z)/H ) 

FZ = ACOS( 1 .0 - 2.0 * ( H(L)-Z(L,I) ) / H(L) ) 
X(LJ) = R(L)/SQP * SQRT( FZ - 0.5*SIN(2*FZ) ) 

33 CONTINUE 



C END OF RADOME LAYER .... 

RETURN 

END 

********************** END op VON KARMAN 

***** ****************************** 

C POWER SERIES SHAPING FUNCTION SUBROUTINE 

SUBROUTINE POWER(L,X,Z,H,R,F,RESOLV,RESPTS ,NLMAX) 

INTEGER IBIG, I, L, SEPRATN, NSHAPE, NLMAX, RESPTS, RESOLV 
REAL X(NLMAX,RESPTS), Z(NLMAX,RESPTS), PI, DELTHETA, DEGTORAD 
REAL THETA,C(5),R(5),H(5),DELZ(5),SHFT,F ACTOR 
REAL ALPHA 1(5), ALPHA(5), B, M, RL.GAMMA, F(NLMAX,7) 

COMPLEX ER.UR 
CHARACTER* 1 YN 

C SET CONSTANTS 
C 

PI=3.14 15927 
DEGTORAD = PI/180. 

C 

C DEFINE START AND END POINT OF EACH RADOME LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RESPTS)=H(L) 

X(L,RESPTS)=0.0 



C DEFINE RADOME LAYER BY OGIVE EQUATIONS .... 
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DELZ(L)=H(L)/RESOLV 



C Prompt for Shaping factor for this power series layer 
C 

WRITER, *) 

WRITE(*,43) 

43 FORMAT(8X,’ The power series shaping function requires a’, 
1/8X, 'shaping factor, F6, to be used in the exponential in', 

1/8X, 'order to determine the radome construction.’, 

1//8X,’ When F6 = 1.0 the radome shape is conical', 

1/8X,' For F6 = 0.75 the radome shape approximates the ', 

1/8X, 'Newtonian contour (minimizes drag). Please enter F6:’) 
READ(*,*) F(L,6) 

DO 33 I=2RESPTS-1 

Z(L,I)=Z(LJ-1)+DELZ(L) 

C EQUATION OF RADOME IS x=R * ( (H-z)/H )**F6 
C 

X(L4) = R(L) * ( (H(L)-Z(L.I)) / H(L) )**F(L,6) 

33 CONTINUE 



C END OF RADOME LAYER .... 

RETURN 

END 

***»****************** END q F powER SERIES 
*********************************** 

C PARABOLIC SHAPING FUNCTION SUBROUTINE 

SUBROUTINE PARACRV(L,X,Z,H,R,F,RESOLV,RESPTS,NLMAX) 

INTEGER IBIG, I, L, SEPRATN, NSHAPE, NLMAX, RESPTS, RESOLV 
REAL X(NLMAX, RESPTS), Z(NLMAXRESPTS), PI, DELTHETA, DEGTORAD 
REAL THETA,C(5),R(5),H(5),DELZ(5),SHFT,F ACTOR 
REAL ALPHA 1(5), ALPHA(5), B, M, RL, GAMMA, F(NLMAX,7) 
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COMPLEX ER,UR 
CHARACTER* 1 YN 

C SET CONSTANTS 
C 

PI=3. 14 15927 
DEGTORAD = PI/180. 

C 

C DEFINE START AND END POINT OF EACH RADOME LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RESPTS)=H(L) 

X(L,RESPTS)=0.0 

C DEFINE RADOME LAYER BY PARABOLIC EQUATION .... 



C Prompt for Shaping factor for this parabolic layer 

C 

WRITE(* ,*) 

WRITE(*,43) 

43 FORMAT(8X,' The parabolic shaping function requires a’, 
l/8X,'shaping factor, F7, to determine the curvature of, 
l/8X,'the parabolic arc for the radome construction.’ 
l//8X,'Please enter F7:') 

READ(*,*)F(L,7) 

C FOR ROUNDED PARABOLIC RADOME LAYER .... 

C ... RADOME IS SCALED TO A HEIGHT = SCALE FUNCTIONS 
DELZ(L)=F(L,7)/RESOLV 
DO 33 I=2,RESPTS-1 
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Z(L,I)=Z(L,I- 1 )+DELZ(L) 



C RECALL EQUATION OF RADOME IS X A 2=F7-Z, X=SQRT(F7-Z) 
X(L,I>SQRT( F(L,7)-Z(L,I) ) 



33 CONTINUE 



C 

C SCALE POINTS TO RADOME DIMENSIONS 
C 

DO 44 I=2,RESPTS-1 
X(L,I)=X(L,I)*R(L)/SQRT( F(L,7) ) 
Z(L,I)=Z(L,I)* H(L)/F(L,7) 

44 CONTINUE 

C END OF ROUNDED RADOME LAYER .... 



C END OF RADOME LAYER .... 

RETURN 

END 

**********************pj^j-) Qp para CURVE 
*********************************** 

C POINTED PARABOLIC SHAPING FUNCTION SUBROUTINE 
SUBROUTINE APEX(L,X,Z,H,R,F,RESOLV,RESPTS,NLMAX) 

INTEGER IBIG, I, L, SEPRATN, NS H APE, NLMAX, RESPTS, RESOLV 
REAL X (NLMAX .RESPTS), Z(NLMAX .RESPTS), PI, DELTHETA, DEGTORAD 
REAL THETA, C(5),R(5),H(5),DELZ(5),SHFT,F ACTOR 
REAL ALPHA 1(5), ALPHA(5), B, M, RL.GAMMA, F(NLMAX,7) 

COMPLEX ER.UR 
CHARACTER* 1 YN 

C SET CONSTANTS 
C 

127 



PI=3. 14 15927 
DEGTORAD = PI/180. 

C 

C DEFINE START AND END POINT OF EACH RADOME LAYER 
C 

Z(L,1)=0.0 

X(L,1)=R(L) 

Z(L,RESPTS )=H(L) 

X(L,RESPTS)=0.0 

C DEFINE RADOME LAYER BY POINTED PARABOLIC EQUATION .... 



C Prompt for Shaping factor for this parabolic layer 
C 

WRITE(*,*) 

WRITE(*,43) 

43 FORMAT(8X,' The APEX parabolic shaping function requires a', 

1/8X, 'shaping factor, F7, to determine the curvature of, 
l/8X,'the parabolic arc for the radome construction and the’, 
l/8X,'percent of the radome that is the pointed apex.’, 

1//8X, 'Please enter F7:') 

READ(*,*)F(L,7) 

C FOR POINTED RADOME LAYER 

C ... RADOME IS SCALED TO A HEIGHT = MODIFIED SCALE FUNCTION ,F7+ 1 
DELZ(L)=( F(L,7)+ 1.0 )/RESOLV 

DO 55 I=2,RESPTS-1 

Z(L,I)=Z(L,I-1)+DELZ(L) 

C 

IF( Z(L,I).LE.( F(L,7) -1 ) ) THEN 
RECALL EQUATION OF CONE IS X A 2=F7-Z, X=SQRT(F7-Z) 
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C 



X(L,I)=SQRT( F(L,7)-Z(L,I) ) 

ELSE 

C STRAIGHT LINE FORMULA FOR POINTED END 
X(LJ)=( ( F(L,7) +1 )-Z(L,I) )/2.0 
ENDIF 



55 CONTINUE 
C 

C SCALE POINTS TO RADOME DIMENSIONS 
C 

DO 66 I=2,RESPTS-1 
X(L,I)=X(L,I)*R(L)/SQRT( F(L,7) ) 
Z(L,I)=Z(L,I)*H(L)/( F(L,7)+1.0 ) 

66 CONTINUE 

C END OF POINTED RADOME LAYER .... 



C END OF RADOME LAYER .... 

RETURN 

END 

********************** END Qp para APEX 
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Appendix C 



CORFLD PROGRAM CODE 



PROGRAM CORFLD 



C Program Created by LT RJ. Vince (Aug-Dec 89) 

C Modifications by MA. Morgan (Nov-Dec 89) 

C Computing core E-field at each point along a 50x50 radome grid 

C in the radome's field coordinate system. 

C 

INTEGER BELL.FF.NPTS ,L J,K,K 1 ,N A.NINC 
INTEGER MSTART.MSTOP.N1 t NMl ( M,MF,MD,Ml,M2,MM 
INTEGER I,ITEJTM,KN,CA,CB»NJNJ s ILP,SELECT,SELCTl 
REAL ALPHA(13),DALPHA( 1 3) ANGINC( 1 3) 

REAL THETA .PHI .PI.DTR .RMIN.RCORE 

REAL XP, YP.ZP.R .COS T.S INT.COS P, S INP.S ,C,C MP.S MP.S A ,S AMX 
REAL R 1 ,R2,RSQ,E M.TH.FI .RPTS.L 1J 1JPLOT(50,50), A(3,3) 

REAL P(40),DP(40),TP(6),PP(6) 

REAL EOTM,EOTE,EI(8),EDTERR(50,50,2),EDTER 1 .EDTER2 
COMPLEX EXP1,E(8),EIX)T(50,50,2),EDOT1,EIX)T2 
COMPLEX EP1 JEP2JET1 .ET2.ER 1 .ER2 J AY,RK,SRE,SRU,KR 
COMPLEX ERAD(50,50,2),ETH(50,50,2),EPHI(50,50,2) 

COMPLEX ER(16),UR(16),COEF(70, 80), CJ(40),DCJ(40),CPLOT 
CHARACTER* 1 YN 
CHARACTER*8 PRTDAT 
CHARACTER* 12 PRTDT1 .PRTDT2.PRTDT3 
CHARACTER*64 GRFLAB.HDR 1 ,HDR2, TITLE 1 ,TITLE2,DUM 
$LARGE 



Setting constants 
BELL=CHAR(7) 

FF = CHAR(12) 

PI = 3.14159265359 
DTR= PI/180. 

JAY= (0.0, 1.0) 

C Set default angle of rotation of planar antenna array away from z axis 

C Set default aspect ratio for plotting 
ASPECT = 1.35 

C Set number of points of resolution across array 
NPTS=50 

RPTS=( FLOAT(NPTS) - 1. )/2. 



C**Initialize arrays 
C Initialize the E-field matrices 
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DO 1 L=1,NPTS 
DO 1 J=1 JVPTS 
DO 1 K=l,2 
ERAD(J r L > K)=(0.,0.) 
ETH(JXJC)=(0.,0.) 
EPHI(J,L,K)=(0.,0.) 

1 CONTINUE 



Prompt for input data filename without extension 



WRITER, 14) 

READ(*,101) PRTDAT 
NCHAR=INDEX(PRTDAT,' ■) -1 
PRTDT1 = PRTDAT(1:NCHAR) // UNT 
OPEN(l 4 JTLE=PRTDT 1 ) 



TITLE l=The datafile used to generate fields is ' 
TITLE 1 =TTTLE 1(1:40) //PRTDT1 

PRTDT2= PRTDAT(1 :NCHAR) // ’.DAT' 
OPEN(l 5 JTLE=PRTDT2) 



Write header from structure data file 

READ(14,101) HDR1 
READ{14,101) HDR2 
READ(14,101) GRFLAB 
WRITEQ 5,101) HDR1 
WRITE(15,101) HDR2 
WRITE(15, 101) GRFLAB 
WRITE!*, 101) HDR1 
WRITE(*,101)HDR2 
WRITER, 101) GRFLAB 



Read interior core Er,Ur 



READO4.507) ER(1),DUM 
READ! 14,507) UR(1),DUM 
SRE=CSQRT(ER(1)) 

SRU=CSQRT(UR(1)) 

KR=SRE*SRU 

WRITER,*) 

WRITE!* ,*)'SRE= ',SRE 
WRITE!*, *)'SRU= ’,SRU 
WRITE!*, *)’KR = ',KR 

WRrrE(15,*)’ SRE=',SRE,'; SRU=',SRU,’; KR = ',KR 



C Read minimum diameter of missile radome core expansion 
C Assuming in Ko*r Wavenumber Units 
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READ(14,503) RMIN.DUM 
RCORE=RMIN 



Read number of incident field angles (.LE. 5) 

READ(14,500) NA.DUM 
NINC=2*NA 



Read incident angles 

IPZ=0 

IMZ=0 

DO 11 1=1 A 

READ(14,5 11) DANGLE 

Checking for +/- Z-Axis Incident Fields 

IF(D ANGLE .EQ.0.0) IPZ=I 
IF(DANGLE.EQ. 180.0) IMZ=I 

DALPHA(I)= DANGLE 

EMRAD uses Poynting vector angle with Z-axis reference 
ANGle of incidence = 180 - DALPHA from EMRAD 

ANGINC(I)= 180- DANGLE 



WRITE(*,*) ' Incoming Inc. Angle #',I 

WRITE(*,*) ' EMRAD Poynting angle of incidence = ’,DALPHA(I) 

WRITE(*,*) ' Radome incident angle = ’,ANGINC(I) 

WRITER,*) 

WRITE(15 *) ’ Incoming Inc. Angle #',1/ = ’.DALPHA(I) 

Set-up array of angles in radian measure 

ALPHA (I)=DTR*DALPHA(I) 

1 1 CONTINUE 
NAI=1 

IF (NA.EQ.l) GO TO 13 
Determine angle of interest 

12 CONTINUE 

WRITE(*,*) ' Select Angle of Interest By Number ' 

WRITER,*) 

READ(*,*)NAI 

IF ((NAI.GE.l).OR.(NAI.LE.NA)) GO TO 13 
WRITE(*,*) ’ That Number is Not Available ... Try Again' 
WRITE(V) 
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GOTO 12 
13 CONTINUE 



C Checking if Selected Field is +/- Z-Axis Incident 
IM1=0 

IF((NAI.EQ.IPZ).OR.(NAI£Q.IMZ)) IM1=1 
WRITE(V) 

WRITE(*,*) ' The angle of consideration is angle #’ ,NAI 
1 = \ANGINC(NAI)/ degrees' 

WRITE(*,*) 

WRITE(15 ,*) ’ The angle of consideration is angle #\NAI 
1 = ’.ANGINCfNAI),’ degrees' 

WRITE(15*) 

C Transformation of coordinates to gimble mounted planar array 
C 

WRITE(*,*)’ The planar array is set by default to lie in the’ 
WRITE(*,*)’ x-y plane, with boresight in the z-direction.' 
WRITER, *) 

WRITE(*,*)’ Do you want to change the array boresight ?’ 
WRITE(*,*)’ ( "Y" or "N” ) ’ 

IFLG=0 
DTH=0.0 
READ(*,101) YN 

IF(YN.EQ.'Y’.OR.YNEQ.'y’) THEN 
WRITE(* *) 

WRITE(*,*) Enter Boresight Theta Angle (in degrees)’ 
READ(V) DTH 
TH = DTH * DTR 
WRITE(* *) 

WRITE(*,*) 'Enter Boresight Phi Angle (in degrees)' 
READ(V) DFI 
FI = DFI * DTR 



IFLG=1 



C Loading Transformation Matrix from Array to Core Coordinates 
C 

CT=COS(TH) 

ST=S1N(TH) 

CP=COS(FI) 

sp=siN(n) 

A(1,1)=CT*CP**2.0+SP**2.0 

A(1,2)=CP*SP*(CT-1.0) 

A(1,3)=ST*CP 

A(2, 1 )=CP* SP*(CT- 1 .0) 

A(2,2)=CT*SP**2.0+CP**2.0 

A(2,3)=ST*SP 

A(3,1)=-ST*CP 
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A(3,2)=-ST*SP 

A(3,3)=CT 

ENDIF 



Get modal information 

READ( 14,505) MSTART.MSTOP.DUM 
WRITEC*,*)’ MSTART=',MSTART,’ MSTOP='MSTOP 
WRITE(15,*) ’ MSTART=',MSTART,' MSTOP=',MSTOP 



Loop modals to read modal coefficients 



WRITE(*,*) ’Reading Coefficients and Constructing’ 

WRITE(*,*) ’Fields for Each Azimuthal "m" Mode* 

WRITEC*,*) 

IF((MSTART.EQ.0).AND.(IM1.EQ.1)) THEN 

WRITE(*,*) Terminating M-Loop at m=l for +/- Z-Axis Incidence’ 

WRITE(*,*) 

MSTOP=l 

ENDIF 

DO 777 MD=MSTART,MSTOP 

READ(14,500) M.DUM 
WRITE(15,*)’ M= ’,M 
EF(M.EQ.MD) GO TO 220 

WRITE(*,*) *Error in Reading "m" - Check Data File’ 

STOP 

220 CONTINUE 
WRITE(*,*) 

WRITEC*,*) ’Azimuthal Mode ”m" = ’,M 
WRITEC*,*) 

IF(M.EQ.0)THEN 
EM = 1.0 
MF = 1 
ELSE 
EM = 2.0 
MF = 0 
ENDIF 

WRITE(*,*) 'Reading Coefficients' 

READ(14,500) NMl.DUM 
DO 665 ID=1,NA 
RE AD( 14,500) I.DUM 
IF(ID.EQ.I) GO TO 221 

WRITE(*,*) 'Inc Angle Index Read Error - Check Data File' 

STOP 

221 CONTINUE 
WRITE(15*)’ 1= ,1 
ITM=I 
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D0 222KN=1,NM1 
CA=KN 
CB=CA+NM1 

READ(14,555) COEF(CA,ITM),COEF(CB,ITM),DUM 
WRITE(15,*) KN,COEF(CA t ITM),COEF(CB,ITM) t DUM 
222 CONTINUE 
ITE=I+NA 
DO 333 KN=1,NM1 
CA=KN 
CB=CA+NM1 

READ(14,555) COEF(CA ,ITE),COEF(CB JTE) ,DUM 
WRITE(15,*) KN,COEF(C A JTE),COEF(CB ,ITE) ,DUM 
333 CONTINUE 
665 CONTINUE 

C Skipping Field Generation for m=0 if +/- Z-Axis Incidence 
IF((M.EQ.0).AND.(IM1.EQ.1)) GO TO 777 



WRITE(*,*) 'Generate Fields Across Array’ 

DX=RCORE/RPTS 

DY=DX 



DO 55 J=1,NPTS 
DO 33 L=1,NPTS 



Determine local planar array coordinates 

J 1 =FLOAT (J- 1 )-RPTS 

Ll=FLOAT(L-l)-RPTS 

XP=J1*DX 

YP=L1*DY 

ZP=0.0 



If Array Tilted Then Transform Coordinates 
From local Planar array coordinates 

to global Core coordinates.... 

IF(IFLG.EQ.l) THEN 
XC= A( 1 , 1 ) * XP+ A( 1 ,2)* YP+ A( 1 ,3)*ZP 
YC=A(2, 1 )* XP+A(2,2)* YP+ A(2,3)*ZP 
ZC=A(3,1)*XP+A(3^)*YP+A(3,3)*ZP 
ELSE 
XC=XP 
YC=YP 
ZC=ZP 
ENDIF 

Determine if R <= RCORE, Else Skip Field Computation 

RSQ=XC*XC+YC*YC+ZC*ZC 

R=SQRT(RSQ) 
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IF (R.GT.RCORE) GO TO 32 



Determine spherical Bessel function scaled radius argument 
RK=KR*R 

Determine Complex Spherical Bessel function and Derivative 

NJN=M+NM 1+MF 

CALL CSBSL(RK,NJN,CJ ,DCJ) 

Compute Core Coordinates and Functions 

THETA=ACOS(ZC/R) 

PHI=ARCTAN2(YC,XC) 

C=COS(THETA) 

S=SIN(THETA) 

CMP=EM*COS(M*PHI) 

SMP=EM*SIN(M*PHI) 

Determine Legendre Polynomial and Derivative 

NLP=NM1+MF 

CALL LPAD(M,NLP,C,P,DP) 

Consider only the incident angle of interest 



I=NAI 

ITM=I 

ITE=ITM+NA 

Initialize additive KN-sum fields to zero 



EP1=(0.,0.) 

EP2=(0.,0.) 

ET1=(0.,0.) 

ET2=(0.,0.) 

ER1=(0.,0.) 

ER2=(0.,0.) 



DO 22 KN=1,NM1 
N=KN+M+MF-1 
NJ=N+1 
NL=KN+MF 
CA=KN 
CB=CA+NM1 

ET 1=ET 1 +COEF(CA,ITM)*M*CJ(NJ)*P(NL) 
1 +COEF(CB,ITM)*Da(NJ)*DP(NL) 

EP 1 =EP 1 +COEF(CA,ITM)*CJ(NJ)* DP(NL) 
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+COEF(CB ,ITM)* M* DCJ(NJ )* P(NL) 
ERl=ERl+COEF(CB,ITM)*N*(N+l)*CJ(NJ)*P(NL) 



ET2=ET2+COEF(CA,ITE)*M*CJ(NJ)*P(NL) 

1 +COEF(CB JTE)* DCJ (NJ)* DP(NL) 

EP2=EP2+COEF(CAJTE)*CJ(NJ)*DP(NL) 

1 +COEF(CBJTE)*M*DCJ(NJ)*P(NL) 

ER2=ER2+COEF(CBJTE)*N*(N+l)*CJ(NJ)*P(NL) 

22 CONTINUE 



ETH(J,L,1)=ETH(J,L 1 1>JAY*SRU*CMP*ET1/R 

EPHI(JJ-,1)=EPHI(JU,1)+JAY*SRU*SMP*EP1/R 

ERAD(J,L,1)=ERAD(J,L ) 1)-JAY*S*CMP*ER1/(SRE*RSQ) 

ETH(J,L,2)=ETH(J,L,2)+SRU*SMP*ET2/R 

EPHI(J,L,2)=EPHI(JU,2)+SRU*CMP*EP2/R 

ERAD(J,L,2)=ERAD(J,L,2)+S*SMP*ER2/(SRE*RSQ) 

32 CONTINUE 

33 CONTINUE 

55 CONTINUE 
C 

C= End of loops across antenna plane ■■ 1 -■■■■■■■= 

777 CONTINUE 
C 

C** End of modal loop for field analysis 
£* ****** ******** 

C Determine magnitude and phase of Computed field in the direction 
C of the theoretical field vector, and the magnitude of the error. 

C 



SINA=SIN(ALPHA(NAI)) 
COS A=COS( ALPHA (NAI» 



C Loop through array to determine field components 
C 

DO 634 J=lJsfPTS 
DO 633 L=1,NPTS 



C 

C 



Determine local planar array coordinates 

Jl=FLOAT(J-l)-RPTS 

Ll=FLOAT(L-l)-RPTS 

XP=J1*DX 

YP=L1*DY 

ZP=0.0 
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C If Array Tilled Then Transform Coordinates 

C From Array to Core .... 

C 

IF(EFLGJEQ.l) THEN 
XC=A(1 ,1)*XP+A(1 ,2)* YP+A(1,3)*ZP 
YC=A(2,1)*XP+A(2,2)*YP+A(2,3)*ZP 
ZC=A(3,1)*XP+A(3,2)*YP+A(3,3)*ZP 
ELSE 
XC=XP 
YC=YP 
ZC=ZP 
ENDIF 

RSQ=XC*XC+YC*YC+ZC*ZC 

R=SQRT(RSQ) 



Determine if R <= RCORE, Else Zero Fields 
IF (R.LE.RCORE) THEN 



Compute Core Coordinates and Functions 

THETA=ACOS(ZC/R) 

PHI=ARCTAN2(YC,XC) 

CT=COS(THETA) 

ST=SIN(THETA) 

CP=COS(PHI) 

SP=SIN(PHI) 

NOTE: R IS GIVEN IN TERMS OF KO*R 

== KO IS INCLUDED IN R, DEFINE KO=l .0 



KO=1.0 



TM and TE Plane waves are of unit magnitude 

EOTE=1.0 

EOTM=1.0 



Complex incident field is exponential(-j*ko*p.r) 
use field direction for incident angle 180 to 90 degrees 

EXP1= CEXP( -JAY*KO*( XC*SINA+ZC*COSA ) ) 

Compute plane wave field unit vectors and components 



E TM = cos(alpha)Ux - sin(alpha)Uz 



C E TM T 
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EI(1) = COSA*CT*CP - SINA*(-ST) 

E(l) = EOTM*EXPl*EI(l) 

C ETMP 

EI(2) = COSA*(-SP) 

E(2) = EOTM*EXPl *EI(2) 

C ETMR 

EI(3) = COSA*ST*CP - SINA*(CT) 

E(3) = EOTM*EXPl*EI(3) 

C E TM TOT 

E(4)=SQRT(CABS(E( 1 ))**2.0+C ABS(E(2))* * 2.0+CABS(E(3))* *2.0) 



E TE = Uy 



ETET 

EI(5) = CT*SP 

E(5) = EOTE*EXPl*EI(5) 

ETEP 
EI(6) = CP 

E(6) = EOTE*EXPl*EI(6) 

ETER 

EI(7) = ST*SP 

E(7) = EOTE*EXPl*EI(7) 

ETETOT 

E(8)=SQRT(CABS(E(5))**2.0+CABS(E(6))**2.0+CABS(E(7))**2.0) 

******** 



Find magnitude and phase of Computed field in the direction 
of the theoretical field vector 

E . ei / Ei = ( eth*ETH+ephi*EPHI+erad*ERAD )/EO*EXPl 

EDOTl=EI(l)*ETH(J^,l)+EI(2)*EPHI(Ja.,l)+EI(3)*ERAD(Jd.,l) 
EDOT (JT-, 1 )=EDOT l/(EOTM*EXP 1 ) 

EDOT2=EI(5)*ETH(J,L,2)+EI(6)*EPHI(J,L > 2)+EI(7)*ERAD(J,L ) 2) 
EDOT(Jd.^)=EDOT2/(EOTE*EXP 1) 

Find rel. magnitude of Computed field perpendicular to the direction 
of the theoretical field vector, the error component in the 
Computed field. 

Rel error = I E - (E.ei) ei I / E0 



= sqrt{ Eth - (E.ei)*ethl**2 + Erad - (E.ei)*eradl**2 
+ Ephi - (E.ei)*ephil**2 )/E0 



Field error magnitude is computed 

EDTER1= CABS(ETH(J,L,1) -EDOT1*EI(1))**2.0 
1 + CABS(EPHI(J,L,1)*EDOT1*EI(2))**2.0 

1 + CABS(ERAD(J^,1)-EDOT1*EI(3))**2.0 
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EDTER2= CABS(ETH(J,L,2) -EDOT2*EI(5))**2.0 
1 + CABS(EPHI(J,L,2)-EDOT2*EI(6))**2.0 

1 + CABS(ERAD(J,L,2)-EDOT2*EI(7))**2.0 

EDTERR(J ,L,2)=S QRT (EDTER 2)/(EOTE) 



ELSE 

C R > RCORE , Set field components to zero 
C 

EDOT(J,L,1)=0.0 
EDOT(J,L,2)=0.0 
EDTERR(J ,L, 1)=0.0 
EDTERR(J,L,2)=0.0 
ENDIF 

633 CONTINUE 

634 CONTINUE 

c********** 



WRITE(* ,*) 'SPHERICAL FIELDS COMPUTED’ 

WRITE(*,*) 

C Optional Listing of Coordinates and E-Fields 
C 

WRITE(*,*) 'Check Values of Coordinates and E-Field ? (Y/N):' 
READ(*,101) YN 

IF((YN.EQ. r Y').OR.(YN.EQ.'y’)) THEN 
757 CONTINUE 

WRITE(*,*) 'Enter Array Triple Index to Check (J,L,I)' 
WRITE(*,*) ’ J= x-point 1,2,..,',NPTS 

WRITE(* ,*) ' L= y-point 1 ,2,../,NPTS 

WRITE(* *) ' 1= 1 for TM, 2 for TE’ 

WRITER,*) 

READ(* ,*) J,L,I 
EPM=CABS(EPHI(J ,L ,1)) 

ETM=CABS(ETH(J,L,I)) 

ERM=CABS(ERAD(JX,I» 

ETOT=SQRT(EPM*EPM+ETM*ETM+ERM*ERM) 



C Determine local planar array coordinates 
C 

Jl=FLOAT(J-l)-RPTS 

Ll=FLOAT(L-l)-RPTS 

XP=J1*DX 

YP=L1*DY 

ZP=0.0 



C If Array Tilted Then Transform Coordinates 
C From Array to Core .... 
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IF(IFLG.EQ.l) THEN 
XC=A(1,1)*XP+A(1,2)*YP+A(1,3)*ZP 
YC=A(2,1)*XP+A(2,2)*YP+A(2,3)*ZP 
ZC=A(3,1)*XP+A(3,2)*YP+A(3,3)*ZP 
ELSE 
XC=XP 
YC=YP 
ZC=ZP 
END IF 

WRITER, *) r XP,YP,ZP:',XP,YP,ZP 
WRITE(*,*) , XC,YC,ZC:',XC,YC,ZC 
WRITE(V) 

WRITE(V) 'IE-RADI =’,ERM 
WRITER, *)'E-THI =’,ETM 
WRITE(*,*) ’E-PHII ='FPM 
WRITE(*,*) ’IE-TOTI =’,ETOT 
WRITE(*,*) 

WRITE(*,*) Look at Another Point ? (Y/N):’ 

READ(*,101) YN 

IF((YN.EQ.T').OR.(YN.EQ.y)) go to 757 
ENDIF 

£*******^******************** ************ ****************************** 

C Field representation selection menu 

p *********************************** 



700 WRITER ,701) 

, ******^* ************************************ ****** ****** ********** 

701 format(////////////5X,’CORFLD OUTPUT GRAPHICS PRESENT ATION'78X, 
l'CORE FIELD REPRESENTATION SELECTION MENUV/5X, 

21*16356 select the type of field representation of interest 
375X, 

4’CORE will display the following field representations:7/8X, 

5*1. TM INCIDENT FIELD, E-THETA',/8X, 

6'2. TM INCIDENT FIELD, E-PHI 78X, 

7’3. TM INCIDENT FIELD, E-RADIAL V8X, 

8'4. TM INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE 7/8X, 

9*5. TE INCIDENT FIELD, E-THETA78X, 

1*6. TE INCIDENT FIELD, E-PHI 78X, 

27. TE INCIDENT FIELD, E-RADIAL 78X, 

3 '8. TE INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE 7/8X, 

4'9. COMPARE COMPUTED FIELD TO INCIDENT FIELD7/8X, 

510. CHANGE ASPECT RATI07/8X, 

6'0. FINISHED WITH THIS ANGLE OF INCIDENCE7/5X, 

6'Indicate your selection entering "10", or ”0" 

7',5X) 

********************************** 



READ(*,*) SELECT 
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C Check response to menu 



IF (SELECT JEQ.l) THEN 

Q ****************************************************************** 

C DISPLAY TM INCIDENT FIELD, E-THETA 

C CHECK FOR DESIRED REPRESENTATION OF FIELD, MAGNITUDE/PHASE 
WRITE(*,715) 

Q ******************* ******************************* **************** 

715 form al(////////////5 X.'CORE FIELD PRESENT ATIONS75X, 
l'You have selected the following field element7/5X, 

2T. TM INCIDENT FIELD, E-THETA7/8X, 

3This field can be plotted in PHASE or MAGNITUDE, V/5X, 

4Please select the type of field representation of interest 
578X, 

6T. PHASE78X, 

72 . MAGNITUDE78X, 

8’Indicate your selection entering "1" or "2" 

97/5X) 

Q ********************************** 

READ(V) SELCT1 

C ASSIGN AS INDICATED BY RESPONSE TO MENU 

DO 712 J=1,NPTS 
DO 711 L=1 J4PTS 
CPLOT=ETH(J,L,l) 

IF (SELCT1.EQ.1) THEN 

PLOT(J,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 

ELSE 

PLOT(J ,L)=C AB S(CPLOT) 

ENDIF 

711 CONTINUE 

712 CONTINUE 
C 

IF (SELCTLEQ.l) THEN 

TITLE2 = ’ TM INCIDENT FIELD, E-THETA PHASE' 

ELSE 

TITLE2 = ' TM INCIDENT FIELD, E-THETA MAGNITUDE' 

ENDIF 



****************************************************************** 



ELSEIF (SELECT.EQ.2) THEN 

Q ****************************************************************** 

C DISPLAY TM INCIDENT FIELD, E-PHI 

C CHECK FOR DESIRED REPRESENTATION OF FIELD, MAGNITUDE/PHASE 
WRITE(*,725) 

/-■ ****************************************************************** 
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725 format(////////////5X,'CORE FIELD PRESENT ATIONS75X, 
l'You have selected the following Field element7/5X, 

2'2. TM INCIDENT FIELD, E-PHIV/8X, 

3This field can be plotted in PHASE or MAGNITUDE, 7/5X, 
4'Please select the type of Field representation of interest 
578X, 

6T. PHASE’,/8X, 

7'2. MAGNITUDE'y8X , 

8'Indicate your selection entering " 1 " or "2" 

97/5X) 

•> ********************************** 



RE ADC ,*) SELCT1 

C ASSIGN AS INDICATED BY RESPONSE TO MENU 



721 

722 
C 



DO 722 J=1,NPTS 
DO 721 L=1,NPTS 
CPLOT=EPHI(J,L,l) 

IF (SELCT1XQ.1) THEN 

PLOT (J ,L)=(ARCT AN2(IMAG(CPLOT) ,REAL(CPLOT)))/DTR 
ELSE 

PLOT(J,L)=CABS(CPLOT) 

ENDIF 

CONTINUE 

CONTINUE 



IF (SELCT1.EQ.1) THEN 
TITLE2 = ’ TM INCIDENT HELD, E-PHI 
ELSE 

TITLE2 = ’ TM INCIDENT HELD, E-PHI 
ENDIF 



PHASE ’ 
MAGNITUDE' 



C 



C 



****************************************************************** 



ELSEIF (SELECT.EQ.3) THEN 

Q ****************************************************************** 

C DISPLAY TM INCIDENT FIELD, E-RADIAL 

C CHECK FOR DESIRED REPRESENTATION OF FIELD, MAGNITUDE/PHASE 
WRITER ,735) 

Q ****************************************************************** 

735 format(////////////5X ,'CORE FIELD PRESENT ATIONS75X, 
l’You have selected the following field element7/5X, 

2'3. TM INCIDENT HELD, E-RADIAL7/8X, 

3 This field can be plotted in PHASE or MAGNITUDE,7/5X, 

4'Please select the type of field representation of interest 
578X, 

6'1. PHASE'yBX, 

7'2. MAGNITUDE' JSX , 

8'Indicate your selection entering " 1 " or "2" 
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97/5X) 

********************************** 



READ(V) SELCT1 

C ASSIGN AS INDICATED BY RESPONSE TO MENU 

DO 732 J=1,NPTS 
DO 731 L=1,NPTS 
CPLOT=ERAD(J,L,l) 

IF (SELCT1.EQ.1) THEN 

PLOT(JU)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 

ELSE 

PLOT(J,L)=CABS(CPLOT) 

ENDIF 

731 CONTINUE 

732 CONTINUE 
C 

IF (SELCT1.EQ.1) THEN 

TITLE2 = ' TM INCIDENT FIELD, E-RADIAL PHASE ’ 

ELSE 

TITLE2 = ' TM INCIDENT FIELD, E-RADIAL MAGNITUDE’ 

ENDIF 

C 

C ****************************************************************** 



ELSEIF (SELECT.EQ.4) THEN 

£ ****************************************************************** 
C DISPLAY TM INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE 

C ASSIGN AS INDICATED BY RESPONSE TO MENU 

DO 742 J=1 J4PTS 
DO 741 L=1,NPTS 
EPM=CABS(EPHI(J JL, 1 )) 

ETM=C ABS (ETH(J JL, 1 )) 

ERM=CABS(ERAD(J,L, 1)) 

PLOT(J,L)=SQRT(EPM*EPM+ETM*ETM+ERM*ERM) 

741 CONTINUE 

742 CONTINUE 

TITLE2 = 1 TM INCIDENT FIELD, E-TOTAL HELD MAGNITUDE ’ 

_ ****************************************************************** 



ELSEIF (SELECT.EQ.5) THEN 

Q ****************************************************************** 

C DISPLAY TE INCIDENT FIELD, E-THETA 

C CHECK FOR DESIRED REPRESENTATION OF FIELD, MAGNITUDE/PHASE 
WRITE(*,755) 

p ****************************************************************** 
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755 format(////////////5X,'CORE FIELD PRESENT ATIONS75X, 
l’You have selected the following field element7/5X, 

2'5. TE INCIDENT FIELD, E-THETAV/8X, 

3 This field can be plotted in PHASE or MAGNITUDE, V/5X, 
41*16386 select the type of field representation of interest 
578X, 

61. PHASE78X, 

72. MAGNITUDE' ^8X , 

81ndicate your selection entering " 1 " or ”2" 

97/5X) 

' ********************************** 



READ(V) SELCT1 

C ASSIGN AS INDICATED BY RESPONSE TO MENU 

DO 752 J=1,NPTS 
DO 751 L=1,NPTS 
CPLOT=ETH(J ,L,2) 

IF (SELCT1.EQ.1) THEN 

PLOT (J ,L)=( ARCT AN2(IMAG(CPLOT),RE AL(CPLOT)))/DTR 
ELSE 

PLOT(J ,L)=C ABS(CPLOT) 

END IF 

751 CONTINUE 

752 CONTINUE 
C 

IF (SELCTlFQ.l)THEN 

TTTLE2 = ' TE INCIDENT FIELD, E-THETA PHASE' 

ELSE 

TITLE2 = ' TE INCIDENT FIELD, E-THETA MAGNITUDE' 

END IF 
C 

c ****************************************************************** 



ELSEIF (SELECT.EQ.6) THEN 

Q ****************************************************************** 

C DISPLAY TE INCIDENT FIELD, E-PHI 

C CHECK FOR DESIRED REPRESENTATION OF FIELD, MAGNITUDE/PHASE 
WRITE(*,765) 

Q ****************************************************************** 

765 format(////////////5X,’CORE FIELD PRESENTATIONS’,^, 
lTou have selected the following field elementV/5X, 

2’6. TE INCIDENT FIELD, E-PHI7/8X, 

3 This field can be plotted in PHASE or MAGNITUDE.7/5X, 

4Please select the type of field representation of interest 
5'JSX, 

61. PHASE78X, 

7’2. MAGN1TUDE78X, 

8’Indicate your selection entering "1" or "2" 
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97/5X) 

********************************** 



READ(* ,*) SELCT1 

C ASSIGN AS INDICATED BY RESPONSE TO MENU 

DO 762 J=1,NPTS 
DO 761 L=1,NPTS 
CPLOT=EPHI(J,L,2) 

IF (SELCT1 JEQ.l) THEN 

PLOT(J,L)=(ARCTAN2(IMAG(CPLOT)JlEAL(CPLOT)))/DTR 

ELSE 

PLOT(J J-)=CABS(CPLOT) 

ENDIF 

761 CONTINUE 

762 CONTINUE 
C 

IF (SELCT1.EQ.1) THEN 

TITLE2 = ' TE INCIDENT FIELD, E-PHI PHASE ’ 

ELSE 

TITLE2 = ’ TE INCIDENT FIELD, E-PHI MAGNITUDE ’ 

ENDIF 

C 

p ****************************************************************** 



ELSEIF (SELECT.EQ.7) THEN 

Q ****************************************************************** 

C DISPLAY TE INCIDENT FIELD, E-RADIAL 

C CHECK FOR DESIRED REPRESENTATION OF HELD, MAGNITUDE/PHASE 
WRITER ,775) 

Q ****************************************************************** 

775 format(////////////5X, , CORE FIELD PRESENTATIONS75X, 

ITou have selected the following field element7/5X, 

27. TE INCIDENT FIELD, E-RADIAL7/8X, 

3This field can be plotted in PHASE or MAGNITUDE, 7/5X, 

4 Please select the type of field representation of interest 
578X, 

61. PHASE78X, 

7'2. MAGNITUDE78X, 

8'Indicate your selection entering "1" or " 2 " 

97/5 X) 

C ********************************** 

READ(*,*) SELCT1 

C ASSIGN AS INDICATED BY RESPONSE TO MENU 



DO 772 J=1,NPTS 
DO 771 L=1 J4PTS 
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CPLOT=ER AD(J JL2) 

IF (SELCT1.EQ.1) THEN 

PLOT(J,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 

ELSE 

PLOT(J ,L)=C ABS (CPLOT) 

ENDIF 

771 CONTINUE 

772 CONTINUE 
C 

IF (SELCTLEQ.l) THEN 

TTTLE2 = ' TE INCIDENT FIELD, E-RADIAL PHASE ' 

ELSE 

TTTLE2 = ’ TE INCIDENT FIELD, E-RADIAL MAGNITUDE' 

ENDIF 

C 

P ****************************************************************** 



ELSEIF (SELECT.EQ.8) THEN 

Q ****************************************************************** 

C DISPLAY TE INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE 

DO 782 J=1,NPTS 
DO 781 L=1 JVPTS 
EPM=C ABS(EPHI(J L-,2)) 

ETM=CABS(ETH(J,L,2)) 

ERM=C AB S(ER AD(J ,L,2)) 

PLOT(J,L)=SQRT(EPM*EPM+ETM*ETM+ERM*ERM) 

781 CONTINUE 

782 CONTINUE 

TITLE2 = ' TE INCIDENT FIELD, E-TOTAL FIELD MAGNITUDE ’ 

****************************************************************** 
ELSEIF (SELECT.EQ.9) THEN 

********************************************************************** 



Incident Field comparison representation selection sub-menu 
************************************************************ 



800 WRITER, 801) 

"« ****************************************************************** 

801 format(////////////5X,'CORE COMPARISON TO IDEAL FIELD7/5X, 
l’CORE FIELD REPRESENTATION SELECTION MENU 27/5X, 

2PIease select the type of field representation of interest 

375X, 

4'CORE will display the following field representations:7/8X, 

51. TM INCIDENT FIELD, Computed Field dot Ideal Field78X, 

6'2. TM INCIDENT FIELD, Scaled Error component 7//8X, 

73. TE INCIDENT FIELD, Computed Field dot Ideal Field',/8X, 

8'4. TE INCIDENT FIELD, Scaled Error component 7//8X, 

9'Indicate your selection entering "1","2","3", or "4" 
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17/5X) 

********************************** 

READ(* ,*) SELCT2 

C Check response to menu 
****** 

IF (SELCT2£Q.l) THEN 

Q ****************************************************************** 

C DISPLAY TM INCIDENT FIELD, Computed Field dot Ideal Field 
C CHECK FOR DESIRED REPRESENTATION OF FIELD, MAGNITUDE/PHASE 
WRITE(*,815) 

Q ****************************************************************** 

815 format(////////////5X,'CORE FIELD PRESENTATIONS75X, 

ITou have selected the following field element7/5X, 

21 . TM INCIDENT FIELD, Computed Field dot Ideal Field7/8X, 

3This field can be plotted in PHASE or MAGNITUDE, 7/5X, 

4’Please select the type of field representation of interest 
578X, 

6T. PHASE78X, 

72. MAGNITUDE’ JSX , 

8'Indicate your selection entering "1" or "2" 

97/5X) 

n ********************************** 



READ(V) SELCT1 

C ASSIGN AS INDICATED BY RESPONSE TO MENU 

DO 812 J=1,NPTS 
DO 811 L=1,NPTS 
CPLOT=EDOT(J,L,l) 

IF (SELCTl£Q.l) THEN 

PLOT(J,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 

ELSE 

PLOT(J,L)=CABS(CPLOT) 

ENDtF 

811 CONTINUE 

812 CONTINUE 
C 

IF (SELCT1.EQ.1) THEN 

TITLE2=TM INCIDENT FIELD, Computed Field dot Ideal Field, PHASE' 
ELSE 

TITLE2=TM INCIDENT FIELD, Computed Field dot Ideal Field, MAGN.’ 
END IF 



****************************************************************** 



ELSEIF (SELCT2.EQ.2) THEN 

****************************************************************** 
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C DISPLAY TM INCIDENT FIELD, Scaled Error component 



DO 822 J=1,NPTS 
DO 821 L=1,NPTS 
PLOT(J,L)=EDTERR(J,L,l) 

821 CONTINUE 

822 CONTINUE 

TITLE2 = TM INCIDENT FIELD, Scaled Error component ’ 

Q ****************************************************************** 

ELSEIF (SELCT2.EQ.3) THEN 

Q ** ************************************************************** 

C DISPLAY TE INCIDENT FIELD, Computed Field dot Ideal Field 
C CHECK FOR DESIRED REPRESENTATION OF FIELD, MAGNITUDE/PHASE 
WRITE(*,815) 

Q ****************************************************************** 

835 format(////////////5X,'CORE FIELD PRESENTATIONS V5X, 
l'You have selected the following field element’ ,//5X, 

2*3. TE INCIDENT FIELD, Computed Field dot Ideal Field7/8X, 

3 This field can be plotted in PHASE or MAGNITUDE,7/5X, 

4Plea.se select the type of field representation of interest 
578X, 

61. PHASE’ ,/8X, 

7'2. MAGNITUDE' J8X , 

8'Indicate your selection entering "1" or ”2" 

97/5X) 

£ ********************************** 

READ(*,*) SELCT1 

C ASSIGN AS INDICATED BY RESPONSE TO MENU 

DO 832 J=1,NPTS 
DO 831 L=1,NPTS 
CPLOT=EDOT(J,L,2) 

IF (SELCT1.EQ.1) THEN 

PLOT(J,L)=(ARCTAN2(IMAG(CPLOT),REAL(CPLOT)))/DTR 

ELSE 

PLOT(J,L)=CABS(CPLOT) 

END IF 

831 CONTINUE 

832 CONTINUE 
C 

IF (SELCT1.EQ.1) THEN 

TITLE2=TE INCIDENT FIELD, Computed Field dot Ideal Field, PHASE’ 

ELSE 

TITLE2=TE INCIDENT FIELD, Computed Field dot Ideal Field, MAGN.' 

END IF 



****************************************************************** 
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ELSEIF (SELCT2.EQ.4) THEN 

Q *************************************************** *************** 

C DISPLAY TE INCIDENT FIELD, Scaled Error component 

DO 842 J=1 J4PTS 
DO 841 L=1 J4PTS 
PLOT(J,L)=EDTERR(J,L,2) 

841 CONTINUE 

842 CONTINUE 

TITLE2 = TE INCIDENT FIELD, Scaled Error component ' 

C ****************************************************************** 

ELSE 

Q ****************************************************************** 

C ( (SELCT21T.O).OR.(SELCT2.GT.4) ) 

GO TO 800 
END IF 



****************** ********************************************** ****** 



ELSEIF (SELECT.EQ.10) THEN 

£ ****************************************************************** 

C ASPECT RATIO 

791 writeC\792) 

792 FORMAT(/5X, '********** ASPECT RATIO **********' J/ 5x, 

TIN ORDER TO AVOID DISTORTION ON THE SCREEN, IT IS NECESSARY75X, 
2TO DEFINE THE NUMBER OF ROWS PER INCH AND THE NUMBER OF75X, 
3'COLUMNS PER INCH. THESE DEFINITIONS CHANGE AS THE75X, 

4 ’RESOLUTION OF THE SCREEN, AS WELL AS THE TYPE OF SCREEN75X, 
5’CHANGE.7/5X, 

C ’THIS PROGRAM WAS WRITTEN IN MODE 16 (EGA) USING75X, 

C 6'A GB-1 VIDEO BOARD AND A NEC MULTISYNC MONITOR. THE75X, 

C 7 'OPTIMUM ASPECT RATIO FOR THIS CONFIGURATION WAS .657//5X, 
8PLEASE INPUT ASPECT RATI0.7/5X) 



read(*,*) aspect 
goto 700 

••* **************** *********************************************** 



ELSEIF (SELECT.EQ.O) THEN 

£ ****************************************************************** 
C FINISHED WITH THIS ANGLE OF INCIDENCE 
goto 999 

£ ****************************************************************** 
ELSE 

£ ****************************************************************** 

C ( (SELECT.LT.O).OR.(SELECT.GT.9) ) 

GO TO 700 
ENDIF 
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Q **********i|^********j|^*************^ *********************************** 



C 



WRITE (*,*) ' CALL PLOTTING ROUTINE' 

CALL PLOT3D(PLOT,50,50,ASPECT > GRFLAB > 

+ TITLE 1 ,TITLE2,HDR 1 ,HDR2,ANGINC(NAI),DTH) 

GO TO 700 



C End of graphic output routine 

999 CONTINUE 
C 

C Program completed 
C 

C I/O FORMAT STATEMENTS 
C 

10 FORMAT (/////////////////7X ,’♦**** WELCOME TO CORFLD *****' J/7X, 

1THIS PROGRAM COMPUTES INTERIOR EM FIELD WITHIN PENETRABLE77X, 
2’BODIES OF REVOLUTION USING THE FIELD COMPONENT INFORMATION',/7X, 
3 'OUTPUT BY THE PROGRAM EMCAD//7X, 

4PLEASE PRESS ANY KEY TO CONTINUE.’/) 

14 FORMAT (//////////////////////7X , 

1THE OUTPUT DATA FILE IS THE FILE CONTAINING THE OUTPUT/7X, 

2 RESULTS OF ALL NUMERICAL CALCULATIONS CONDUCTED BY EMCAD.’./7X, 
3 THE FORMAT FOR THIS INPUT IS FILENAME ONLY. NO EXTENSION IS’/7X, 

4 REQUIRED OR DESIRED. EMCAD AUTOMATICALLY APPENDS AN ’/7X, 

5 'EXTENSION OF .DMT TO YOUR FILENAME.7/7X, 

6PLEASE ENTER THE FILENAME OF THE OUTPUT DATA FILE.7/7X) 

23 FORMAT(/7X,ENTER INC FLD ANGLE (DEG) FOR # ’.I3//7X) 

25 FORMAT (//////7X, PLEAS E PRESS ANY KEY TO CONTINUE.'////) 

100 FORMAT(/ '******************** CORE OUTPUT DATA **************** 

l*****»**»**'j 

101 FORMAT(A) 

102 FORMAT(I5) 

103 FORMAT(4(E14.6)) 

104 FORMAT(' *** PROGRAM ABORTED BECAUSE **«") 

107 F0RMAT(//7X,I5,2(2X,’(',1PE11.3,2X,1PE1 1.3,')’)) 

109 FORMAT(//7X,'INCIDENT FIELD ANGLES7/9X,’N’,10X,’ALPHA(N)') 

1 10 FORMAT(//7X,I3,7Xp5.0,' DEG 1 ) 

111 FORMAT(/7X, 'COMPLEX Er(n) AND Ur(n) 1 ) 

1 12 FORMAT (JflX .PNTER INTERIOR DATA FILE (D:FILENAME.EXTENSION): ’) 

1 15 FORMAT(//7X, r ENTER CAPTION OR LABEL ') 

124 FORMAT(//7X, 'INDEXING PROGRAM THROUGH VALUES OF M') 

125 FORMAT(7X,'M-LOOP .... M = ’,15) 

135 FORMAT(7X,PX M-LOOP’) 

1 36 FORMAT(//7X, ’********* CORE ANALYSIS COMPLETED **************•) 

204 FORMAT(//7X,' ************** M = ’,15,' **************'//) 

300 FORMAT(7X,I5,5E12.3) 



uuu * uuuuuuuuu uu 



500 FORMAT(3X,I5,A64) 

503 FORMAT(3X,E14.6,A64) 

505 FORMAT(3X,I5,2XJ5A64) 

507 FORMAT(3X,lPEl 1.3,2X,1PE1 1.3,2X,A64) 

510 FORMAT(3X,F5.0.A64) 

511 FORMAT(3X,F5.0) 

555FORMAT(3X,E12.63X,E12.6,4X,E12.6,3X,E12.6,3X,A64) 

566 FORMAT(3X,A10^ 12.6,3X^ 12.6) 

567 FORMAT(3X,A 103 12.6) 

576 FORMAT(3X,A5,I5A2312.6,3X312.6) 

577 FORMAT(3X,A5,E12.6) 

578 FORMAT(3X312.6) 



99 CONTINUE 

Program termination 

WRITER, 136) 

STOP 

END 

****** ******************************************************* 



SUBROUTINE CSBSL (Z.Nl.CJ.DCJ) 

COMPLEX Z32.F1 32 .FSTO, ALPHA, CJ(* ),DCJ(*) 

This routine computes Riccati spherical bessel functions and 
derivatives for complex Z using downward recurrence as in 
Miller's algorithm, but with highly accurate starting values 
using truncated Tatlor's series. 

Returns CJ(n)=J[n-l] (kr) for n=l...Nl 
DCJ(n)=J'[n-l] (kr) 



DEFINING CONSTANTS 

ZM2=(CABS(Z))**2 
N MX=N 1 +ZM2 
Z2=(Z*Z)/2.0 
Dl=2.*NMX+3. 

D2=Dl*(2.*NMX+5.) 

D3=D2*(2.*NMX+7.) 

D4=D3*(2.*NMX+9.) 

C Using Truncated Taylor's series for relative upper starting 
C values of J(NMX) and J(NMX+1) with NMX .GT. ABS(Z*Z). 
C 

F1=1.-Z2/D1+Z2*Z2/(2.*D2)-Z2*Z2*Z2/(6.*D3) 
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F2=Z*(1 7D1 -Z2/D2+Z2*Z2/(2.*D3)-Z2*Z2*Z2/(6.*D4)) 
Performing downward recurrence with highest two orders in CJ(N1). 



M=NMX 
1 1 FSTO=Fl 
F 1=(2.*M+ 1 ,)*F1/Z-F2 
F2=FST0 

IF (C ABS(F1 ).LT. 1 .0E20) GO TO 1 
Fl=Fl*1.0E-20 
F2=F2* 1 .0E-20 
FSTO=FSTO* 1 .0E-20 
1 CONTINUE 
M=M-1 

IF (M+2.GT.N1) GO TO 1 1 
CJ(N1)=F2 
CJ(N1-1)=F1 
N0=Nl-2 



Continuing downward recurrence to CJ(1) 



DO 22 K=1,N0 
J=N1-K-1 

22 CJ(J)=(2.*J+l.)*CJ(J+l)/Z-a(J+2) 



Normalizing entire array wrt actual J0(Z). 

ALPHA=CSIN(Z)/CJ(1) 

DO 33 K=l r Nl 
33 CJ (K)= ALPHA *CJ (K) 



Computing derivatives 

DCJ ( 1 )=CJ ( 1 )/Z-CJ (2) 

DO 44 K=2Js t 1 
44 DCJ (K)=CJ (K- 1 )-(K- 1 )* CJ (K)/Z 
RETURN 
END 

*********************************** 

SUBROUTINE LPAD (MJM1,C,P,DP) 

C 

C This routine computes Normalized Legendre polynomials and 
C derivatives. 

C 

C For m >= 1; P(n) = P(m,n+m-l)/sin(THETA) 

C DP(n) = dP(m,n+m-l) / dTHETA for n=l,2,... 

C For m = 0; P(n) = P(0,n)/sin(THETA) 

C DP(n) = dP(0,n) / dTHETA for n=l,2,... 

C 



REAL P(*),DP(*) 
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KSP=N1-1 

Q=1.0-C*C 

S=SQRT(Q) 

IF(M.GT.O) GO TO 1 

P(l)=1.0 

P(2)=C 

DP(1)=0.0 

DP(2)=-S 

Q=S 

GO TO 3 

1 IF (M.EQ.l) GO TO 5 
P(l)=-((-S)**(M-l)) 

DP(1)=-M*C*(-S)**(M-1) 

GO TO 7 

5 P(l)=-1.0 
DP(1)=-M*C 
7 DO 2 K=1,M 
P(1)=P(1)*(2*K-1) 

2 DP( 1 )=DP( 1 )* (2* K- 1) 

P(2)=P(1)*C*(2*M+1) 

DP(2)=(C* DP( 1 )-Q*P( 1 ))*(2* M+ 1 ) 

3 DO 4 K=2,KSP 
N=K+M-1 
AK=M-N-1 

P(K+ 1 )=((M+N)* P(K- 1 )-(2* N+ 1 )*C* P(K))/AK 

4 DP(K+1)=((M+N)*DP(K-1)+(2*N+1)*(Q*P(K)-C*DP(K)))/AK 
RETURN 

END 

FUNCTION ARCTAN2(Y,X) 



Routine returns ATAN(Y/X) for all cases 

REAL ARCTAN2,Y,X r PI 
PI=3. 1415927 
IF(Y.NE.0.0) GO TO 11 
IF(X.EQ.0.0) GO TO 33 
IF(X.GT.0.0) ARCT AN 2=0.0 
IF(X.LT.0.0) ARCTAN2=PI 
RETURN 
11 IF(X.NE.0.0) GO TO 22 
IF(Y.GT.0.0) ARCTAN2=0.5*PI 
IF(Y.LT.0.0) ARCTAN2=1.5*PI 
RETURN 

22 ARCTAN2=ATAN2(Y,X) 

RETURN 
33 ARCTAN2=0.0 
RETURN 
END 

c * * * % ** * * * * ** * * * * *** * ** * * * * * * ** * ** ** *** ** **** * * *** ** * ** ** * *** * * * * * 
C PLOT3D Plotting subroutine 
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C REVISED: 22 NOV 1989 
C 

C This subroutine automatically scales and plots data as 3-D surfaces 
C 

C All subroutines beginning with 'Q' 

C are from the Grafmatic plotting library 

C 

SUBROUTINE PLOT3D(Z,M,N,ASPECT,GRFLAB, 

+ TITLE 1 ,TTTLE2 ,HDR 1 ,HDR2,ALPHA,G AMMA) 

REAL X(50,50),Y(50,50),Z(50,50),XMIN,XMAX,YMIN,YMAX,ZMAX,ZMIN 
REAL XST,XFIN,YST,YFTN,ZST r ZFIN r ZFLOOR,ZAVG 
REAL Pffl,THETA t ALPHA,GAMMA,STHETA,IUl > M2 1 N2,R 
INTEGER M.N 

CHARACTER* 1 YES ,GR AFPRNT 

CHARACTER*5 ASTRJNG,GSTRING,TSTRING,PSTRING,MIN,MAX,AVG 
CHARACTER* 12 name, label 

CHARACTER*64 TITLE 1, TITLE2, TITLE3, TITLE4, GRFLAB, HDR1, HDR2 

CHARACTER * 64 HCOPY 

DIMENSION IWORK 1 (640),IWORK2(640) 

C SETTING CONSTANTS 
C 

BELL=CHAR(7) 

FF = CHAR(12) 



C Set constant for degree to radian conversion 
C 

PI = 3.14159265359 
DTR= PI/180. 

C Calculate # of points on array 
C 

NARRAY=PI*( (M/2)* *2) 

M2=FLOAT(M-1)/2.0 

N2=FLOAT(N-1)/2.0 

C IWORK1 and IWORK2 should have a dimension equal to the number of 
C pixels (or points on a pen plotter) in the x direction. 

C 

NPTS=640 

C NPTS is the size of the IWORK1 and IWORK2 arrays 
C 

MODEP=16 

C Plot will be in graphics mode 16 
C 

C XYSET is a routine called to set up the initial x,y values (see below) 

C 

CALL XYSET(X,Y,M,N) 
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XST = X(l,l) 
YST = Y(l,l) 
ZST = Z(l,l) 

XFIN = XST 
YFIN = YST 
ZFEN = ZST 



ZMIN = 1000 
ZMAX = -1000 
ZAVG = 0.0 

PHI = 45.0 
THETA = 75.0 



DO 10 I = 1,M 
DO 15 J = 1,N 

IF (X(IJ) .LT. XST) XST = X(U) 

IF (X(U) .GT. XFIN) XFIN = X(U) 

IF (Y(U) .LT. YST) YST = Y(IJ) 

IF (Y(U) .GT. YFIN) YFIN = Y(U) 

IF (Z(U) .LT. ZST) ZST = Z(U) 

IF (Z(U) .GT. ZFIN) ZFIN = Z(U) 



C Check if point is on the antenna array 
C 

I1=FLOAT(I)-1.0 

J1=FLOAT(J)-1.0 

R= SQRT( ((Il-M2)/M2)**2.0 + ((Jl-N2)/N2)**2.0 ) 
IF (R.LT.1.0) THEN 
IF (Z(U) .LT. ZMIN) ZMIN = Z(IJ) 

IF (Z(U) .GT. ZMAX) ZMAX = Z(U) 

ZAVG = ZAVG + Z(U)/(FLOAT(N ARRAY)) 
ENDIF 

15 CONTINUE 
10 CONTINUE 
XFIN = XFIN + .2 
YFIN = YFIN + .2 



C The following parameters set the axes values, etc. (see documentation 
C for Q3DXAX, Q3DYAX, Q3DZAX) 

C 

XBEG=XST 

YBEG=YST 

XEND=XFIN 

YEND=YFIN 

XSTEP = .2 
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YSTEP = .2 

IF ((XFIN-XST).GE. 0.9) XSTEP = .5 
IF ((YFTN-YST).GE. 0.9) YSTEP = .5 



ZDIF=ZFTN-ZST 



IF (ZDIF.LE. 0.3) THEN 
ZSTEP = 0.05 

ZFIN = FLOAT(INT(10*ZFIN))/10+2*ZSTEP 
ZST = FLOAT(INT(10*ZST))/10-ZSTEP 
ZMINOR = 1 

ELSEIF (ZDIFXE. 1 .0) THEN 
ZSTEP = 0.25 

ZFIN = FLOAT(INT(ZFIN+0.5))+2*ZSTEP 
ZST = FLOAT(INT(ZST-0.5))-ZSTEP 
ZMINOR = 1 

ELSEIF (ZDIF.LE. 2.0) THEN 
ZSTEP = 0.5 

ZFIN = FLOAT(INT(ZFIN+0.5))+2*ZSTEP 
ZST = FLOAT(INT(ZST-0.5))-ZSTEP 
ZMINOR = 1 

ELSEIF (ZDIFXE. 5.0) THEN 
ZSTEP = 1.0 

ZFIN = FLOAT(INT(ZFIN+0.5))+2*ZSTEP 
ZST = FLOAT(INT(ZST-0.5))-ZSTEP 
ZMINOR = 1 

ELSEIF (ZDIFXE. 10.0) THEN 
ZSTEP = 2.0 

ZFIN = FLOAT(INT(ZFIN+0.5))+2*ZSTEP 
ZST = FLOAT(INT(ZST-0.5))-ZSTEP 
ZMINOR = 1 
ELSE 

ZSTEP = 30.0 

ZFIN = INT (ZFIN +ZSTEP) 

ZST = INT(ZST-15.0) 

ZMINOR = 1 
ENDIF 
ZFLOOR=0.0 



Clear the screen and enter input data 

CALL QSMODE(2) 

C 3D-stick is a "wire-frame" figure with hidden lines removed. 

C 3D-fill is similar, except that the interior regions of each 
C "panel" making up the surface if filled in with color "klrin" 

C and the boundary of each panel is drawn with color "klredg". 

C NOTE: The two routines use different hidden line removal algorithm. 
C Which you use is a matter of taste. However the equivalent 

C of Q3DFIL is not available on the pen plotter since it is 

C based on an overpainting scheme which isn't practical on a plotter. 
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c 

C The angles phi and theta specify the view direction as shown in the 
C appropriate figure in your GRAFMATIC/PLOTMATIC manual. 
C The parameter icross is one if you wish cross-grid lines in both the 
C x and y directions or zero for grid lines only parallel to 

C the x axis (as may be approrpiate in creating a spectral plot). 

C 

C 

16 WRITE(*,*) ' Enter 0 for 3D-stick or 1 for 3D-fill option' 

C READ(V) IFIL 

IFIL=0 
WRITE(* *) 

WRITE(*,*) ' Default view angles are: phi=-45.,theta=70.' 

IF (IFIL.EQ.1) THEN 

KLRIN is the color inside the panel 

KLRIN=2 

KLREDG is the panel border color 

KLREDG=1 

ELSE 

ICROSS is 1 for cross-grids in both directions 

ICROSS=l 
END IF 

Continuation point to loop through desired angles of observation 

17 CONTINUE 

CALL QSMODE(MODEP) 

Modify THETA to reflect YFTN/ZFTN 

STHETA = ATAN( (YFIN/ZFTN)*TAN(THETA*DTR) )/DTR 



Q3DROT changes the view perspective as described in the documentation 

CALL Q3DROT(X,Y,Z,M,N ( PHI,STHETA) 

C 

CALL Q3DWIND(XST,XFIN,YST,YFIN,ZST,ZFIN,XMIN,XMAX,YMIN,YMAX) 
IOPT=0 

IF ( (XMAX-XMIN).NE.0.0) THEN 
YOVERX = (YMAX-YMIN) / (XMAX-XMIN) 

ELSE 

YOVERX = 999 
ENDIF 
XORG=0.0 
YORG=0.0 
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CALL QPLOT(100,550,30,270,XMIN,XMAX,YMIN,YMAX, 
1 XORG.Y ORG.IOPT, Y OVERX, ASPECT) 



C The next 3 routines create the x, y, and z axes 
C If Q3DFIL is used, then draw the axes AFTER the 3D plot is drawn 
C so that the 3D plot does not overwrite the x/y axes. 

C 

IF (IFIL.EQ.l) GO TO 66 
C 

CALL Q3DXAX(XST,XFIN,XSTEP,9,1,2,YBEG,YEND,ZFLOOR,1.0) 
CALL Q3DYAX(YST,YFIN,YSTEP,9,1,2,XBEG,XEND,ZFLOOR,1.0) 
CALLQ3DZAX(ZST,ZFIN,ZSTEP,9,-1,2,XBEG,YBEG,1.0) 

66 CONTINUE 



Finally draw the 3D plot, using the appropriate option. 

IF(IFIL.EQ.l) CALL Q3DFIL(X,Y,M,N,KLRIN,KLREDG) 

IF(IFIL.NE.l) CALL Q3DSTk(X,Y,M,N,IWORKl,rWORK2,NPTS,ICROSS) 

Draw the axes if using Q3DFIL routine 

IF(IFIL.NE.l) GO TO 67 

CALL Q3DXAX(XST,XFIN,XSTEP,9,1,2,YBEG,YEND,ZFLOOR,1.0) 
CALLQ3DYAX(YST,YFIN,YSTEP,9,1,2,XBEG,XENDZFLOOR,1.0) 

CALL Q3DZAX(ZST,ZFIN,ZSTEP,9,- 1 ,2,XBEG, YBEG, 1 .0) 

67 CONTINUE 

Q ****************************************************************** 

C WRITE TEXT ON SCREEN 
NCHAR = 64 
KOLOR = 10 
ICOL = 10 
IROW = 24 

CALL QPTXT (NCHAR, TITLE 1 ,KOLOR,lCOL,IROW) 

£ ****************************************************************** 
C WRITE TEXT ON SCREEN 
NCHAR = 64 
KOLOR = 10 
ICOL= 10 
IROW = 23 

CALL QPTXT(NCHAR,HDR1 , KOLOR, ICOL,IROW) 

£ ****************************************************************** 
C WRITE TEXT ON SCREEN 
NCHAR = 64 
KOLOR = 10 
ICOL = 10 
IROW = 22 

CALL QPTXT(NCHAR,HDR2,KOLOR,ICOL,IROW) 

Q ****************************************************************** 

C WRITE TEXT ON SCREEN 
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NCHAR = 64 
KOLOR = 10 
ICOL = 10 
IROW = 21 

CALL QPTXT(NCHAR,GRFLAB,KOLOR,ICOL,IROW) 

Q ****************************************************************** 

C WRITE TEXT ON SCREEN 
NCHAR = 64 
KOLOR = 10 
ICOL= 15 

mow = 20 

CALL QPTXT (NCHAR ,TITLE2,KOLOR JCOL.IROW) 

Q ****************************************************************** 

C MAGNITUDES INTO A STRING 
C 

CALL NUMSTR(ZMIN,MIN) 

CALL NUMSTR(ZMAX.MAX) 

CALL NUMSTR(Z A VG ,AVG) 

TITLE3=’min. Z = '//MIN//; max. Z = '//MAX 
+ //; average Z = ’//AVG 
C WRITE TEXT ON SCREEN 
C 

NCHAR = 64 
KOLOR = 10 
ICOL = 17 
mow = 3 

CALL QPTXT(NCHAR,TITLE3,KOLOR,ICOL,mOW) 

( 2 ****************************************************************** 

C ANGLES INTO A STRING 
C 

CALL NUMSTR(ALPHA,ASTRING) 

CALL NUMSTR(GAMMA,GSTRING) 

CALL NUMSTR(THETA,TSTRING) 

CALL NUMSTR(PHI,PSTRING) 

TITLE4='lnc. A = '//ASTRING// , ;Plane G = ’//GSTRING 
+ //; View: phi='//PSTRING//' lheta='//TSTRING 
C WRITE TEXT ON SCREEN 
C 

NCHAR = 64 
KOLOR = 10 
ICOL = 10 

mow = 2 

CALL QPTXT(NCHAR,TITLE4,KOLORJCOL,mOW) 

Q ****************************************************************** 

C Hardcopy Query 
NCHAR = 30 
KOLOR = 12 
ICOL = 1 
IROW = 1 

HCOPY=’Hardcopy — > Enter P or p’ 

CALL QPTXT(NCHAR,HCOPY, KOLOR, ICOL.IROW) 

CALL QCMOV(32,l) 
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CALL QINKEY (EXTEND, IKEY) 

NCHAR = 40 
KOLOR = 1 
ICOL = 1 
IROW = 1 
HCOPY=' 

CALL QPTXT(NCHAR JiCOPY.KOLOR.ICOLJROW) 

C If P* or 'p' then prtsc 
C 

IF( IKEY.EQ.80 .OR. IKEY.EQ.l 12 ) THEN 
CALL QPSCRN 
END IF 

Q ****************************************************************** 

CALL Q3DINV(X,Y,Z,M ,N) 

CALL QSMODE(2) 

WRITE(* ,*)'ENTER NEW ANGLES PHI & THETA (DEG) OR (9,9) WHEN DONE' 
RE AD(5,*)PHI, THETA 

IF ( (PHI.NE.9.).AND.(THETA.NE.9.) ) GOTO 17 
999 CONTINUE 

100 FORMAT(Al) 

520 FORMAT(3X,I5) 

578 FORMAT(3X,E12.6) 

101 FORMAT(A) 

RETURN 

END 



***************************** 

Q ************************************************************* 

C 

SUBROUTINE XYSET(X,Y,M,N) 

C 

REAL X(M,N),Y(M,N),I1,J1,N2,M2 
C 

C define x,y values over the m by n grid 
C 

M2=FLOAT(M-1)/2.0 
N2=FLOAT(N-1)/2.0 
DO 111 1=1, M 
DO 121 J=1,N 
I1=FLOAT(I)-1.0 
J1=FLOAT(J)-1.0 
X(IJ)= (Il-M2)/M2 
Y(IJ)= (Jl-N2)/N2 
121 CONTINUE 
111 CONTINUE 
RETURN 
END 

Q* *************** ***************** 

C 

C SUBROUTINE NUMSTR 
C 
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ACCEPTS NUMBER AND RETURNS CHAR*5 WITH SIGN AND THREE DIGIT STRING 

SUBROUTINE NUMSTR(ANUM.STR) 

REAL ANUM.VNUM 

INTEGER NUM.HUN.TEN.ONE.ORDER 

CHARACTER’S STR 

ORDER = 0 



Get sign information 

IF (ANUM.GE.0.0) THEN 
STR(1:1)=V 
ELSE 

STR(1:1)=’-' 

END IF 

Strip sign values 

VNUM=ABS(ANUM) 

IF (VNUM .LT. 1 .0) THEN 
ORDER=0 

VNUM=VNUM* 1000.0 
ELSEIF (VNUM .LT .10.0) THEN 
ORDER=l 

VNUM= VNUM* 100.0 
ELSEIF (VNUM .LT. 100.0) THEN 
ORDER=2 
VNUM= VNUM* 10.0 
ELSE 
ORDER=3 
ENDIF 

Use integer part of shifted VNUM 



NUM=VNUM 

IF (NUM.GE.1000) THEN NUM=999 
HUN=INT(NUM/100.) 
TEN=(NUM-100*HUN)/10 
ONE=(NUM-100*HUN-10*TEN) 

PLACE DECIMAL 

IF (ORDER .EQ. 0) THEN 
STR(2:2)='.' 

STR(3:3)=CHAR(48+HUN) 

STR(4:4)=CHAR(48+TEN) 

STR(5:5)=CHAR(48+ONE) 

ELSEIF (ORDER .EQ. 1)THEN 
STR(2:2)=CHAR(48+HUN) 
STR(3:3)='.' 
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STR(4:4)=CHAR(48+TEN) 
STR(5:5)=CHAR(48+ONE) 
ELSEIF (ORDER .EQ. 2) THEN 
STR(2:2)=CHAR(48+HUN) 
STR(3:3)=CHAR(48+TEN) 
STR(4:4)=’.’ 

STR(5:5)=CHAR(48+ONE) 

ELSE 

STR(2:2)=CHAR(48+HUN) 

STR(3:3)=CHAR(48+TEN) 

STR(4:4)=CHAR(48+ONE) 

STR(5:5)=7 

ENDIF 



RETURN 

END 
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SOFTWARE SOURCES 



Microsoft FORTRAN 

16011 NE 36th Way 
Box 97017 

Redmond, WA 98073-9717 

West Coast Consultants CURVE DIGITIZER 
4202 Genesee Avenue, Suite 309 
San Diego, CA 92117 
(619)565-1266 

Microcompatibles, Inc. GRAFMATIC 
301 Prelude Drive 
Silver Spring, MD 20901 
(301)593-5151 

Prof. M.A. Morgan 

Code 62Mw 

Naval Postgraduate School 
Monterey, CA 93943 
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